← Back to branch summary

~maddevelopers/mg5amcnlo/WWW5_caching

« back to all changes in this revision

Viewing changes to users/mardelcourt/PROC_141512/PROC_141512/SubProcesses/projection.f

  • Committer: John Doe
  • Date: 2013-03-25 20:27:02 UTC
  • Revision ID: john.doe@gmail.com-20130325202702-5sk3t1r8h33ca4p4
first clean version

Show diffs side-by-side

added added

removed removed

Lines of Context:
users/mardelcourt/PROC_141512/PROC_141512/SubProcesses/projection.f
 
1
c++++++++++++++++++++++
 
2
c  CONTENT
 
3
c++++++++++++++++++++++
 
4
 
 
5
c      real*8  function PROJ_ETA(P,NHEL,IC)
 
6
c      real*8 function PROJ_psi(P,NHEL,IC)
 
7
c      real*8 function PROJ_PSI_REL(P,NHEL,IC)
 
8
c      real*8 function PROJ_CHI(P,NHEL,IC,J_qn)
 
9
c      real*8 function PROJ_H(P,NHEL,IC)
 
10
c      real*8 function PROJ_PSILEPT(P,NHEL,IC)
 
11
 
 
12
c      subroutine spin_projection(P,NHEL,IC,TEMP)
 
13
c      subroutine spin_singlet_proj(P,NHEL,IC,TEMP)
 
14
c      subroutine pseudoscalar(P1,P2,LAMBDA1,LAMBDA2,MC,PS)
 
15
c      subroutine pseudoscalar0(P1,P2,LAMBDA1,LAMBDA2,M1,M2,PS)
 
16
c      subroutine vector0(P1,P2,LAMBDA1,LAMBDA2,M1,M2,jio)
 
17
c      subroutine VECTOR(P1,P2,LAMBDA1,LAMBDA2,MQ,jio)
 
18
c      subroutine vectorBc(P1,P2,LAMBDA1,LAMBDA2,M1,M2,jio)
 
19
c      subroutine leptonic_current(k1,k2,l1,l2,mq, jio)
 
20
c      subroutine genmom_lept(P,x,k1,k2,three_k1)
 
21
c      SUBROUTINE polariz(P,EPS)
 
22
c      subroutine get_helicity(P,eps)
 
23
c      subroutine get_helicity2(P,eps)
 
24
 
 
25
 
 
26
 
 
27
      real*8  function PROJ_ETA(P,NHEL,IC)
 
28
C
 
29
C******************************************************************************
 
30
C     THIS FUNCTION PROJECTS THE AMPLITUDE ONTO A CC~(3S1) STATE               * 
 
31
C     AND SQUARE IT                                                           *
 
32
C     THE ROUTINE IS UNIVERSAL                                                *
 
33
C******************************************************************************
 
34
C
 
35
      implicit none
 
36
      include 'nexternal.inc'
 
37
      include 'coupl.inc'
 
38
      include 'color_data.inc'
 
39
C
 
40
C     PARAMETERS
 
41
C
 
42
      INTEGER  NEXT_AMP                    ! number of  particles
 
43
      PARAMETER (NEXT_AMP=NEXTERNAL+1)     ! # of external particles before the projection
 
44
C
 
45
C     ARGUMENTS
 
46
C
 
47
      double precision  P(0:3,NEXTERNAL)
 
48
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
 
49
C
 
50
C     LOCAL VARIABLES
 
51
C
 
52
C         1. for the projection
 
53
      INTEGER NHEL_AMP(next_amp),IC_AMP(next_amp), L1,L2,J,K,n
 
54
      DOUBLE COMPLEX PS
 
55
      DOUBLE PRECISION P_AMP(0:3,next_amp),
 
56
     &                  P1(0:3),P2(0:3)  ! MOMENTA OF c AND c~ 
 
57
     & , weight1,weight2
 
58
C
 
59
      COMPLEX*16 JAMP(NCOLOR)
 
60
      DOUBLE COMPLEX PROJAMP(NCOLOR)
 
61
C
 
62
C         2.  to  square the amplitude
 
63
      COMPLEX*16 ZTEMP
 
64
C
 
65
c--masses
 
66
      DOUBLE PRECISION       qmass(2)
 
67
      COMMON/to_qmass/qmass      
 
68
C
 
69
C----
 
70
C  Begin code
 
71
C----
 
72
C     SPECIFYING THE MOMENTA OF C,C~ + THE TOTAL MOMENTUM
 
73
C
 
74
 
 
75
      weight1=qmass(1)/(qmass(1)+qmass(2))
 
76
      weight2=qmass(2)/(qmass(1)+qmass(2))
 
77
 
 
78
      DO L1=0,3
 
79
      DO L2=1, next_amp-2
 
80
      P_AMP(L1,L2)=P(L1,L2)
 
81
      ENDDO
 
82
      P_AMP(L1,next_amp-1)=P(L1,nexternal)*weight1
 
83
      P_AMP(L1,next_amp)=P(L1,nexternal)*weight2
 
84
      P1(L1)=P(L1,nexternal)*weight1
 
85
      P2(L1)=P(L1,nexternal)*weight2
 
86
      ENDDO
 
87
C
 
88
C     RECORD THE POLARIZATIONS
 
89
C
 
90
      DO J=1,next_amp-2
 
91
      NHEL_AMP(J)=NHEL(J)
 
92
      IC_AMP(J)=IC(J)
 
93
      ENDDO
 
94
      IC_AMP(next_amp-1)=1
 
95
      IC_AMP(next_amp)=1
 
96
C
 
97
C************************************* 
 
98
C     START THE SUM   
 
99
C************************************* 
 
100
C
 
101
C     INITIALIZATION
 
102
C
 
103
      DO K=1,NCOLOR
 
104
       projamp(K)=(0.0D0,0.0D0)
 
105
      ENDDO
 
106
 
 
107
      DO L1=-1,1,2  ! SUM OVER HELICITIES
 
108
      DO L2=-1,1,2
 
109
      NHEL_AMP(next_amp-1)=L1
 
110
      NHEL_AMP(next_amp)=L2 
 
111
c      write(*,*) 'P1',P1
 
112
c      write(*,*) 'P2',P2
 
113
c      write(*,*) 'L1',L1
 
114
c      write(*,*) 'L2',L2
 
115
      CALL pseudoscalar0(P1,P2,L1,L2,qmass(1),
 
116
     & qmass(2),PS)
 
117
c      write(*,*) 'PS',PS
 
118
      CALL matrix(P_AMP,NHEL_AMP,IC_AMP,JAMP)
 
119
c      write(*,*) 'Jamp',Jamp
 
120
C 
 
121
      DO K=1, NCOLOR    ! INDEX OF JAMP
 
122
          projamp(K)=projamp(K)+JAMP(K)*PS  ! THE COLOR FACTOR IS INCLUDED IN SQAMP
 
123
      ENDDO
 
124
C  
 
125
      ENDDO
 
126
      ENDDO
 
127
C       
 
128
C************************************ 
 
129
C    SQUARE THE AMPLITUDE
 
130
C***********************************
 
131
C
 
132
      PROJ_ETA = 0.D0
 
133
      DO n = 1, NCOLOR
 
134
          ZTEMP = (0.D0,0.D0)
 
135
          DO J = 1, NCOLOR
 
136
              ZTEMP = ZTEMP + CF(J,n)*PROJAMP(J)
 
137
          ENDDO
 
138
          PROJ_ETA =PROJ_ETA
 
139
     &         +dble(ZTEMP*DCONJG(PROJAMP(n))/DENOM(n))
 
140
      ENDDO
 
141
      RETURN 
 
142
      END 
 
143
C
 
144
      real*8 function PROJ_psi(P,NHEL,IC)
 
145
C
 
146
C******************************************************************************
 
147
C     THIS FUNCTION PROJECT THE AMPLITUDE ONTO A QQ~(3S1) STATE               * 
 
148
C     AND SQUARE IT                                                           *
 
149
C     THE ROUTINE IS UNIVERSAL                                                *
 
150
C******************************************************************************
 
151
C
 
152
      implicit none
 
153
      include 'nexternal.inc'
 
154
      include 'coupl.inc'
 
155
      include 'color_data.inc'
 
156
C
 
157
C     PARAMETERS
 
158
C
 
159
      INTEGER  NEXT_AMP                   
 
160
      PARAMETER (NEXT_AMP=NEXTERNAL+1)     ! # of external particles before the projection
 
161
C
 
162
C     ARGUMENTS
 
163
C
 
164
      REAL*8 P(0:3,NEXTERNAL)
 
165
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
 
166
C
 
167
C     LOCAL VARIABLES
 
168
C
 
169
C         1. for the projection
 
170
      INTEGER NHEL_AMP(next_amp),IC_AMP(next_amp), L1,L2,J,K
 
171
      DOUBLE COMPLEX JIO(6),TEMP(4,NCOLOR)
 
172
      DOUBLE PRECISION P_AMP(0:3,next_amp),
 
173
     &                  P1(0:3),P2(0:3),PTOT(0:3)   ! MOMENTA OF c AND c~ AND TOTAL MOMENTUM
 
174
      DOUBLE complex EPS(0:3,-1:1)
 
175
C
 
176
      COMPLEX*16 JAMP(NCOLOR)
 
177
      DOUBLE COMPLEX PROJAMP(NCOLOR)
 
178
C
 
179
C         2.  to  square the amplitude
 
180
      COMPLEX*16 ZTEMP
 
181
C
 
182
c--masses
 
183
      DOUBLE PRECISION       qmass(2)
 
184
      COMMON/to_qmass/qmass      
 
185
      DOUBLE PRECISION weight1,weight2
 
186
C
 
187
C----
 
188
C  Begin code
 
189
C----
 
190
C     SPECIFYING THE MOMENTA OF C,C~ + THE TOTAL MOMENTUM
 
191
C
 
192
 
 
193
      weight1=qmass(1)/(qmass(1)+qmass(2))
 
194
      weight2=qmass(2)/(qmass(1)+qmass(2))
 
195
 
 
196
      DO L1=0,3
 
197
      DO L2=1, next_amp-2
 
198
      P_AMP(L1,L2)=P(L1,L2)
 
199
      ENDDO
 
200
      P_AMP(L1,next_amp-1)=P(L1,nexternal)*weight1
 
201
      P_AMP(L1,next_amp)=P(L1,nexternal)*weight2
 
202
      P1(L1)=P(L1,nexternal)*weight1
 
203
      P2(L1)=P(L1,nexternal)*weight2
 
204
      PTOT(L1)=P1(L1)+P2(L1)
 
205
      ENDDO
 
206
 
 
207
C
 
208
C     RECORD THE POLARIZATIONS
 
209
C
 
210
      DO J=1,next_amp-2
 
211
      NHEL_AMP(J)=NHEL(J)
 
212
      IC_AMP(J)=IC(J)
 
213
      ENDDO
 
214
      IC_AMP(next_amp-1)=1
 
215
      IC_AMP(next_amp)=1
 
216
C
 
217
C************************************* 
 
218
C     START THE SUM   
 
219
C************************************* 
 
220
C
 
221
C     INITIALIZATION
 
222
C
 
223
      DO J=1,4,1
 
224
      DO K=1,NCOLOR
 
225
       TEMP(J,K)=(0.0D0,0.0D0)
 
226
      ENDDO
 
227
      ENDDO
 
228
 
 
229
      DO L1=-1,1,2  ! SUM OVER HELICITIES
 
230
      DO L2=-1,1,2
 
231
      NHEL_AMP(next_amp-1)=L1
 
232
      NHEL_AMP(next_amp)=L2  
 
233
      CALL VECTOR0(P1,P2,L1,L2,qmass(1),qmass(2),JIO)
 
234
 
 
235
 
 
236
      CALL matrix(P_AMP,NHEL_AMP,IC_AMP,JAMP)
 
237
C 
 
238
C 
 
239
      DO K=1, NCOLOR    ! INDEX OF JAMP
 
240
          DO J=1,4,1   ! PAY ATTENTION TO THE RANGE : J IN [1,4]
 
241
          TEMP(J,K)=TEMP(J,K)+JAMP(K)*JIO(J)  ! THE COLOR FACTOR IS INCLUDED IN SQAMP
 
242
          ENDDO
 
243
      ENDDO
 
244
C 
 
245
      ENDDO
 
246
      ENDDO
 
247
C       
 
248
c       CALL polariz(PTOT,EPS)
 
249
       CALL get_helicity2(PTOT,EPS,(qmass(1)+qmass(2)))
 
250
 
 
251
C
 
252
      DO K=1,NCOLOR
 
253
      PROJAMP(K)=TEMP(1,K)*EPS(0,NHEL(nexternal))-
 
254
     &  TEMP(2,K)*EPS(1,NHEL(nexternal))-
 
255
     & TEMP(3,K)*EPS(2,NHEL(nexternal))- 
 
256
     & TEMP(4,K)*EPS(3,NHEL(nexternal))
 
257
      ENDDO
 
258
C     
 
259
 
 
260
C************************************ 
 
261
C    SQUARE THE AMPLITUDE
 
262
C***********************************
 
263
C
 
264
      PROJ_PSI = 0.D0
 
265
      DO I = 1, NCOLOR
 
266
          ZTEMP = (0.D0,0.D0)
 
267
          DO J = 1, NCOLOR
 
268
              ZTEMP = ZTEMP + CF(J,I)*PROJAMP(J)
 
269
          ENDDO
 
270
          PROJ_PSI =PROJ_PSI
 
271
     &         +DBLE(ZTEMP*DCONJG(PROJAMP(I)))/DENOM(I)
 
272
      ENDDO
 
273
      RETURN 
 
274
      END 
 
275
C
 
276
      real*8 function PROJ_PSI_REL(P,NHEL,IC)
 
277
C
 
278
C******************************************************************************
 
279
C     THIS FUNCTION PROJECT THE AMPLITUDE ONTO A CC~(3S1) STATE               * 
 
280
C     AND SQUARE IT                                                           *
 
281
C     THE ROUTINE IS UNIVERSAL                                                *
 
282
C******************************************************************************
 
283
C
 
284
      include 'nexternal.inc'
 
285
      include 'coupl.inc'
 
286
      include 'color_data.inc'
 
287
C
 
288
C     PARAMETERS
 
289
C
 
290
      INTEGER  NEXT_AMP                    ! number of  particles
 
291
      PARAMETER (NEXT_AMP=NEXTERNAL+1)     ! # of external particles before the projection
 
292
C
 
293
C     ARGUMENTS
 
294
C
 
295
      REAL*8 P(0:3,NEXTERNAL)
 
296
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
 
297
C
 
298
C     LOCAL VARIABLES
 
299
C
 
300
C         1. for the projection
 
301
      INTEGER NHEL_AMP(next_amp),IC_AMP(next_amp), L1,L2,J,K
 
302
      DOUBLE COMPLEX JIO(6),TEMP(4,NCOLOR)
 
303
      DOUBLE COMPLEX JIOSTAR(6),TEMPSTAR(4,NCOLOR)
 
304
      DOUBLE PRECISION P_AMP(0:3,next_amp),
 
305
     &                  P1(0:3),P2(0:3),PTOT(0:3)   ! MOMENTA OF c AND c~ AND TOTAL MOMENTUM
 
306
      DOUBLE PRECISION P_AMPSTAR(0:3,next_amp),
 
307
     &                  P1STAR(0:3),P2STAR(0:3)
 
308
      DOUBLE PRECISION EPS(0:3,-1:1)
 
309
C
 
310
      COMPLEX*16 JAMP(NCOLOR)
 
311
      DOUBLE COMPLEX PROJAMP(NCOLOR)
 
312
      COMPLEX*16 JAMPSTAR(NCOLOR)
 
313
      DOUBLE COMPLEX PROJAMPSTAR(NCOLOR)
 
314
C
 
315
C         2.  to  square the amplitude
 
316
      COMPLEX*16 ZTEMP
 
317
      complex*16 proj_matrix_temp
 
318
C
 
319
c--masses
 
320
      DOUBLE PRECISION       qmass(2)
 
321
      COMMON/to_qmass/qmass      
 
322
C
 
323
C     global
 
324
C
 
325
      double precision prel1(0:3),prel2(0:3)
 
326
      double precision BOUND, RHOSQ          ! BOUND IS THE UPPER BOUND FOR q_{REL}^2
 
327
      common /to_rel_mom/prel1,prel2,BOUND, RHOSQ
 
328
C
 
329
C----
 
330
C  Begin code
 
331
C----
 
332
C     SPECIFYING THE MOMENTA OF C,C~ + THE TOTAL MOMENTUM
 
333
C
 
334
      DO L1=0,3
 
335
      DO L2=1, next_amp-2
 
336
      P_AMP(L1,L2)=P(L1,L2)
 
337
      P_AMPSTAR(L1,L2)=P(L1,L2)
 
338
      ENDDO
 
339
C                 in the amplitude
 
340
      P_AMP(L1,next_amp-1)=P(L1,nexternal)/2+prel1(L1)
 
341
      P_AMP(L1,next_amp)=P(L1,nexternal)/2-prel1(L1)
 
342
      P1(L1)=P(L1,nexternal)/2+prel1(L1)
 
343
      P2(L1)=P(L1,nexternal)/2-prel1(L1)
 
344
C                 in the compl. conj. amplitude 
 
345
      P_AMPSTAR(L1,next_amp-1)=P(L1,nexternal)/2+prel2(L1)
 
346
      P_AMPSTAR(L1,next_amp)=P(L1,nexternal)/2-prel2(L1)
 
347
      P1STAR(L1)=P(L1,nexternal)/2+prel2(L1)
 
348
      P2STAR(L1)=P(L1,nexternal)/2-prel2(L1)
 
349
c
 
350
      PTOT(L1)=P(L1,nexternal)
 
351
      ENDDO
 
352
C
 
353
C     RECORD THE POLARIZATIONS
 
354
C
 
355
      DO J=1,next_amp-2
 
356
      NHEL_AMP(J)=NHEL(J)
 
357
      IC_AMP(J)=IC(J)
 
358
      ENDDO
 
359
      IC_AMP(next_amp-1)=1
 
360
      IC_AMP(next_amp)=1
 
361
C
 
362
C************************************* 
 
363
C     START THE SUM   
 
364
C************************************* 
 
365
 
 
366
C
 
367
C     INITIALIZATION
 
368
C
 
369
      DO J=1,4,1
 
370
      DO K=1,NCOLOR
 
371
       TEMP(J,K)=(0.0D0,0.0D0)
 
372
       TEMPSTAR(J,K)=(0.0D0,0.0D0)
 
373
      ENDDO
 
374
      ENDDO
 
375
 
 
376
      DO L1=-1,1,2  ! SUM OVER HELICITIES
 
377
      DO L2=-1,1,2
 
378
      NHEL_AMP(next_amp-1)=L1
 
379
      NHEL_AMP(next_amp)=L2  
 
380
      CALL VECTOR(P1,P2,L1,L2,qmass(1),JIO)              ! for amp.
 
381
      CALL VECTOR(P1STAR,P2STAR,L1,L2,qmass(1),JIOSTAR)  ! for amp*
 
382
 
 
383
      CALL matrix(P_AMP,NHEL_AMP,IC_AMP,JAMP)         ! for amp.
 
384
      CALL matrix(P_AMPSTAR,NHEL_AMP,IC_AMP,JAMPSTAR) ! for amp*
 
385
C 
 
386
c       
 
387
      DO K=1, NCOLOR    ! INDEX OF JAMP
 
388
          DO J=1,4,1   ! PAY ATTENTION TO THE RANGE : J IN [1,4]
 
389
C                  THE COLOR FACTOR IS INCLUDED IN SQAMP
 
390
          TEMP(J,K)=TEMP(J,K)+JAMP(K)*JIO(J)  ! for amp
 
391
          TEMPSTAR(J,K)=TEMPSTAR(J,K)+JAMPSTAR(K)*JIOSTAR(J)  ! for amp*
 
392
          ENDDO
 
393
      ENDDO
 
394
C 
 
395
C  
 
396
      ENDDO
 
397
      ENDDO
 
398
C       
 
399
      CALL polariz(PTOT,EPS)
 
400
C
 
401
      DO K=1,NCOLOR
 
402
c                          for amp.
 
403
      PROJAMP(K)=TEMP(1,K)*EPS(0,NHEL(nexternal))-
 
404
     &  TEMP(2,K)*EPS(1,NHEL(nexternal))-
 
405
     & TEMP(3,K)*EPS(2,NHEL(nexternal))- 
 
406
     & TEMP(4,K)*EPS(3,NHEL(nexternal))
 
407
C                           for amp*
 
408
      PROJAMPSTAR(K)=TEMPSTAR(1,K)*EPS(0,NHEL(nexternal))-
 
409
     &  TEMPSTAR(2,K)*EPS(1,NHEL(nexternal))-
 
410
     & TEMPSTAR(3,K)*EPS(2,NHEL(nexternal))-
 
411
     & TEMPSTAR(4,K)*EPS(3,NHEL(nexternal))
 
412
      ENDDO
 
413
C     
 
414
 
 
415
C************************************ 
 
416
C    SQUARE THE AMPLITUDE
 
417
C***********************************
 
418
C
 
419
      PROJ_MATRIX_TEMP = (0.0D0,0.0D0)
 
420
      DO I = 1, NCOLOR
 
421
          ZTEMP = (0.D0,0.D0)
 
422
          DO J = 1, NCOLOR
 
423
              ZTEMP = ZTEMP + CF(J,I)*PROJAMP(J)
 
424
          ENDDO
 
425
          PROJ_MATRIX_TEMP =PROJ_MATRIX_TEMP
 
426
     &         +ZTEMP*DCONJG(PROJAMPSTAR(I))/DENOM(I)
 
427
      ENDDO
 
428
      
 
429
      PROJ_PSI_REL = DBLE(proj_matrix_temp)
 
430
      RETURN 
 
431
      END 
 
432
C
 
433
      real*8 function PROJ_CHI(PP,NHEL,IC,J_qn)
 
434
C
 
435
C******************************************************************************
 
436
C     THIS FUNCTION PROJECT THE AMPLITUDE ONTO A P wave CC~ STATE             * 
 
437
C     AND SQUARE IT                                                           *
 
438
C     THE ROUTINE IS UNIVERSAL.                                               *     
 
439
C     MOMENTA MUST BE GIVEN IN THE QUARKONIUM REST FRAME                      *
 
440
C     UNDER WORK         !!                                                   *
 
441
C******************************************************************************
 
442
C
 
443
      implicit none
 
444
      include 'nexternal.inc'
 
445
      include 'coupl.inc'
 
446
      include 'color_data.inc'
 
447
C
 
448
C     ARGUMENTS
 
449
C
 
450
      REAL*8 PP(0:3,NEXTERNAL)
 
451
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL),J_qn
 
452
C
 
453
C     LOCAL VARIABLES
 
454
C
 
455
C         1. for the projection
 
456
      DOUBLE PRECISION Pboost(0:3)  !,Pboostnd(0:3)
 
457
      Double precision P(0:3,NEXTERNAL),Z(0:3)
 
458
      double complex eps(1:3,-1:1)
 
459
c      double precision Ponium(0:3),Poniumnd(0:3)
 
460
C
 
461
      DOUBLE COMPLEX TEMP(4,NCOLOR),
 
462
     & TEMP1(4,NCOLOR),TEMP2(4,NCOLOR),TEMP3(4,NCOLOR)
 
463
C
 
464
      DOUBLE COMPLEX AMP(3,3,NCOLOR)
 
465
      DOUBLE COMPLEX hel_AMP(NCOLOR)
 
466
      DOUBLE COMPLEX hel_AMPJm1(NCOLOR,3),spin_mat(NCOLOR,3,-1:1)
 
467
C
 
468
C         2.  to  square the amplitude
 
469
      double precision sqm2,sqm6,sqm3
 
470
      COMPLEX*16 ZTEMP
 
471
      integer j,k,l,ia,ib,mu
 
472
C
 
473
c--masses
 
474
      DOUBLE PRECISION       qmass(2)
 
475
      COMMON/to_qmass/qmass
 
476
C
 
477
C     global
 
478
C
 
479
      double precision prel(0:3)
 
480
      common /to_pWAVE/prel
 
481
c
 
482
      double precision pnd(0:3,nexternal),delta
 
483
 
 
484
      double complex rho11,rho00,rho1m1,rhom1m1,rhom11
 
485
      double precision plam,pnu
 
486
      common /to_pol_param/rho11,rho00,rho1m1,rhom1m1,rhom11,plam,pnu
 
487
 
 
488
C
 
489
C----
 
490
C  Begin code
 
491
C----
 
492
      delta=0.00001d0
 
493
      
 
494
c
 
495
c   go to the quarkonium rest frame
 
496
c
 
497
 
 
498
 
 
499
       pboost(0)=pp(0,nexternal)
 
500
       do k=1,3
 
501
       pboost(k)=-pp(k,nexternal)
 
502
       enddo
 
503
       do j=1,nexternal
 
504
       call boostx(PP(0,j),pboost,P(0,j))
 
505
       do mu=0,3
 
506
       Pnd(mu,j)=P(mu,j)
 
507
       enddo
 
508
       enddo
 
509
       Pnd(0,nexternal)=dsqrt(P(0,nexternal)**2+delta**2)
 
510
c
 
511
c    define the tensor projected amplitude with UPPER indices
 
512
c
 
513
c    PAY ATTENTION : d/dq_i = - d/d_q^j
 
514
c
 
515
c      a) no relative momentum
 
516
c
 
517
       do k=0,3
 
518
       prel(K)=0d0
 
519
       enddo
 
520
c
 
521
       prel(1)=delta/1000000d0
 
522
       call spin_projection(P,NHEL,IC,TEMP)
 
523
c
 
524
c       a) component q_1
 
525
c
 
526
       prel(1)=delta
 
527
 
 
528
       call spin_projection(Pnd,NHEL,IC,TEMP1)
 
529
c
 
530
c       b) component q_2
 
531
c
 
532
       prel(1)=0d0
 
533
       prel(2)=delta
 
534
 
 
535
       call spin_projection(Pnd,NHEL,IC,TEMP2)
 
536
c
 
537
c       c) component q_3
 
538
c
 
539
       prel(2)=0d0
 
540
       prel(3)=delta
 
541
 
 
542
       call spin_projection(Pnd,NHEL,IC,TEMP3)
 
543
c
 
544
c      computing derivatives
 
545
c
 
546
       do k=1,NCOLOR
 
547
         do j=1,3
 
548
c       amp(0,j,k)=(0d0, 0d0)
 
549
       amp(1,j,k)=-(temp1(j+1,k)-temp(j+1,k))/dcmplx(delta)
 
550
       amp(2,j,k)=-(temp2(j+1,k)-temp(j+1,k))/dcmplx(delta)
 
551
       amp(3,j,k)=-(temp3(j+1,k)-temp(j+1,k))/dcmplx(delta)
 
552
         enddo
 
553
       enddo
 
554
 
 
555
C************************************ 
 
556
C    Select helicity state
 
557
C***********************************
 
558
C
 
559
 
 
560
      call get_helicity(PP(0,nexternal),eps)
 
561
 
 
562
       sqm6=1d0/dsqrt(6d0)
 
563
       sqm2=1d0/dsqrt(2d0)
 
564
       sqm3=1d0/dsqrt(3d0)
 
565
 
 
566
      if(j_qn.eq.0) then
 
567
 
 
568
      do L=1,Ncolor
 
569
      hel_amp(L)=(0d0,0d0)
 
570
 
 
571
          do Ia=1,3
 
572
          do IB=1,3
 
573
           hel_amp(l)=hel_amp(l)+ amp(Ia,Ib,L)*(eps(Ia,1)*eps(Ib,-1)
 
574
     &+eps(Ia,-1)*eps(Ib,1) -eps(Ia,0)*eps(Ib,0))*sqm3 
 
575
          enddo
 
576
          enddo
 
577
      enddo !L
 
578
 
 
579
 
 
580
      elseif(j_qn.eq.1) then
 
581
 
 
582
      do L=1,Ncolor
 
583
      hel_amp(L)=(0d0,0d0)
 
584
 
 
585
        if (nhel(nexternal).eq.1 ) then
 
586
          do Ia=1,3
 
587
          do IB=1,3
 
588
           hel_amp(l)=hel_amp(l)+ 
 
589
     & amp(Ia,Ib,L)*(eps(Ia,1)*eps(Ib,0)-eps(Ia,0)*eps(Ib,1))*sqm2 
 
590
          enddo
 
591
          enddo
 
592
 
 
593
        elseif (nhel(nexternal).eq.0 ) then
 
594
          do Ia=1,3
 
595
          do IB=1,3
 
596
           hel_amp(l)=hel_amp(l)+ 
 
597
     & amp(Ia,Ib,L)*(eps(Ia,1)*eps(Ib,-1)-eps(Ia,-1)*eps(Ib,1))*sqm2 
 
598
          enddo
 
599
          enddo
 
600
 
 
601
        elseif (nhel(nexternal).eq.-1 ) then
 
602
          do Ia=1,3
 
603
          do IB=1,3
 
604
           hel_amp(l)=hel_amp(l)+
 
605
     & amp(Ia,Ib,L)*(eps(Ia,0)*eps(Ib,-1)-eps(Ia,-1)*eps(Ib,0))*sqm2
 
606
          enddo
 
607
          enddo
 
608
        endif
 
609
 
 
610
      enddo !L
 
611
 
 
612
      elseif(j_qn.eq.2) then
 
613
 
 
614
      do L=1,Ncolor
 
615
      hel_amp(L)=(0d0,0d0)
 
616
 
 
617
        if (nhel(nexternal).eq.2 ) then
 
618
          do Ia=1,3
 
619
          do IB=1,3
 
620
           hel_amp(l)=hel_amp(l)+ 
 
621
     & amp(Ia,Ib,L)*(eps(Ia,1)*eps(Ib,1)) 
 
622
          enddo
 
623
          enddo
 
624
 
 
625
        elseif (nhel(nexternal).eq.1 ) then
 
626
          do Ia=1,3
 
627
          do IB=1,3
 
628
           hel_amp(l)=hel_amp(l)+ 
 
629
     & amp(Ia,Ib,L)*(eps(Ia,1)*eps(Ib,0)+eps(Ia,0)*eps(Ib,1))*sqm2 
 
630
          enddo
 
631
          enddo
 
632
 
 
633
        elseif (nhel(nexternal).eq.0 ) then
 
634
          do Ia=1,3
 
635
          do IB=1,3
 
636
           hel_amp(l)=hel_amp(l)+
 
637
     & amp(Ia,Ib,L)*(eps(Ia,1)*eps(Ib,-1)+eps(Ia,-1)*eps(Ib,1)+
 
638
     & 2d0*eps(Ia,0)*eps(Ib,0))*sqm6
 
639
          enddo
 
640
          enddo
 
641
 
 
642
        elseif (nhel(nexternal).eq.-1 ) then
 
643
          do Ia=1,3
 
644
          do IB=1,3
 
645
           hel_amp(l)=hel_amp(l)+
 
646
     & amp(Ia,Ib,L)*(eps(Ia,0)*eps(Ib,-1)+eps(Ia,-1)*eps(Ib,0))*sqm2
 
647
          enddo
 
648
          enddo
 
649
 
 
650
        elseif (nhel(nexternal).eq.-2 ) then
 
651
          do Ia=1,3
 
652
          do IB=1,3
 
653
           hel_amp(l)=hel_amp(l)+
 
654
     & amp(Ia,Ib,L)*eps(Ia,-1)*eps(Ib,-1)
 
655
          enddo
 
656
          enddo
 
657
        endif
 
658
     
 
659
c      write(*,*) 'hel_amp for hel,col',nhel(nexternal),L,hel_amp(l)
 
660
      enddo !L
 
661
 
 
662
 
 
663
      elseif(j_qn.eq.-1) then
 
664
c     here is a special option: 
 
665
c       * I average over the angular momentum degrees of freedom
 
666
c       * I keep track of the spin of the QQ~ pair
 
667
c     This is used for polarization studies involving transitions 3PJ[8]> J/psi
 
668
      do L=1,Ncolor
 
669
      do Ia=1,3
 
670
      hel_ampJm1(L,Ia)=(0d0,0d0)
 
671
c
 
672
          do IB=1,3
 
673
           hel_ampJm1(l,Ia)=hel_ampJm1(l,Ia)+
 
674
     & amp(Ia,Ib,L)*eps(Ib,nhel(nexternal))
 
675
          enddo
 
676
 
 
677
      enddo !angular momentum
 
678
      enddo !color
 
679
 
 
680
      ZTEMP = (0.D0,0.D0)
 
681
      do Ia=1,3 !angular momentum
 
682
      DO L = 1, NCOLOR
 
683
          DO K = 1, NCOLOR
 
684
          ZTEMP = ZTEMP+CF(K,L)*hel_ampJm1(K,Ia)*dconjg(hel_ampJm1(L,Ia))/DENOM(L)
 
685
          ENDDO
 
686
      ENDDO
 
687
      ENDDO
 
688
      PROJ_CHI = DBLE(ztemp)
 
689
      return
 
690
c
 
691
      elseif(j_qn.eq.-7) then
 
692
c     here is another special option:
 
693
c       * I average over the angular momentum degrees of freedom
 
694
c       * I compute the spin density matrix
 
695
c     This is used for polarization studies involving transitions 3PJ[8]> J/psi
 
696
      do L=1,Ncolor
 
697
      do Ia=1,3
 
698
        hel_ampJm1(L,Ia)=(0d0,0d0)
 
699
        spin_mat(l,Ia,-1)=(0d0,0d0)
 
700
        spin_mat(l,Ia,1)=(0d0,0d0)
 
701
        spin_mat(l,Ia,0)=(0d0,0d0)
 
702
c
 
703
          do IB=1,3
 
704
 
 
705
           spin_mat(l,Ia,-1)=spin_mat(l,Ia,-1)+
 
706
     & amp(Ia,Ib,L)*eps(Ib,-1)
 
707
           spin_mat(l,Ia,1)=spin_mat(l,Ia,1)+
 
708
     & amp(Ia,Ib,L)*eps(Ib,1)
 
709
           spin_mat(l,Ia,0)=spin_mat(l,Ia,0)+
 
710
     & amp(Ia,Ib,L)*eps(Ib,0)
 
711
          enddo
 
712
          enddo
 
713
      enddo !L
 
714
 
 
715
      ZTEMP = (0.D0,0.D0)
 
716
      do Ia=1,3
 
717
      do Ib=1,3
 
718
      DO L = 1, NCOLOR
 
719
          DO K = 1, NCOLOR
 
720
          ZTEMP = ZTEMP+CF(K,L)*amp(Ia,Ib,K)*dconjg(amp(Ia,Ib,L))/DENOM(L)
 
721
          ENDDO
 
722
      ENDDO
 
723
      ENDDO
 
724
      ENDDO
 
725
      PROJ_CHI = DBLE(ztemp)
 
726
 
 
727
      do Ia=1,3
 
728
      DO L = 1, NCOLOR
 
729
          DO K = 1, NCOLOR
 
730
          rho11 = rho11+CF(K,L)*spin_mat(K,Ia,1)*dconjg(spin_mat(L,Ia,1))/DENOM(L)
 
731
          rho00 = rho00+CF(K,L)*spin_mat(K,Ia,0)*dconjg(spin_mat(L,Ia,0))/DENOM(L)
 
732
          rho1m1 = rho1m1+CF(K,L)*spin_mat(K,Ia,1)*dconjg(spin_mat(L,Ia,-1))/DENOM(L)
 
733
          rhom11 = rhom11+CF(K,L)*spin_mat(K,Ia,-1)*dconjg(spin_mat(L,Ia,1))/DENOM(L)
 
734
          rhom1m1 = rhom1m1+CF(K,L)*spin_mat(K,Ia,-1)*dconjg(spin_mat(L,Ia,-1))/DENOM(L)
 
735
          ENDDO
 
736
      ENDDO
 
737
      ENDDO
 
738
 
 
739
      return
 
740
 
 
741
 
 
742
      endif ! end if J_qn= ...
 
743
 
 
744
C************************************ 
 
745
C     sum over colored amplitudes
 
746
C***********************************
 
747
 
 
748
      ZTEMP = (0.D0,0.D0)
 
749
      DO L = 1, NCOLOR
 
750
          DO K = 1, NCOLOR
 
751
          ZTEMP = ZTEMP+CF(K,L)*hel_amp(K)*dconjg(hel_amp(L))/DENOM(L)
 
752
          ENDDO
 
753
c          PROJ_chi =proj_chi+ZTEMP*DCONJG(JAMP(L))/DENOM(L)
 
754
      ENDDO
 
755
      PROJ_CHI = DBLE(ztemp)
 
756
c      write(*,*) 'proj_chi',proj_chi
 
757
c      pause 
 
758
 
 
759
c
 
760
c
 
761
c      Here I directly sum over the chi_c polarization (not used anymore)
 
762
c
 
763
c      ztemp=(0d0,0d0)
 
764
c      do Ia=2,4
 
765
c      do IB=2,4
 
766
c      DO K=1,ncolor
 
767
c      DO L=1,ncolor
 
768
c      if (j_qn.eq.0) then
 
769
c      ztemp=ztemp+amp(Ia,Ia,K)*DCONJG(amp(Ib,Ib,L))/3d0
 
770
c     &  *CF(K,L)/DENOM(L)
 
771
c       elseif (j_qn.eq.1) then
 
772
c      ztemp=ztemp+(amp(Ia,Ib,K)*DCONJG(amp(Ia,Ib,L))-
 
773
c     &     amp(Ia,Ib,K)*DCONJG(amp(Ib,Ia,L)))
 
774
c     &  *CF(K,L)/DENOM(L)/2d0
 
775
c      elseif (j_qn.eq.2) then
 
776
c      ztemp=ztemp+((amp(Ia,Ib,K)*DCONJG(amp(Ia,Ib,L))+
 
777
c     &     amp(Ia,Ib,K)*DCONJG(amp(Ib,Ia,L)))/2d0-
 
778
c     &   amp(Ia,Ia,K)*DCONJG(amp(Ib,Ib,L))/3d0)
 
779
c     &  *CF(K,L)/DENOM(L)
 
780
c      endif
 
781
c      enddo
 
782
c      enddo
 
783
c      enddo
 
784
c      enddo
 
785
c
 
786
c      PROJ_CHI = DBLE(ztemp)
 
787
 
 
788
      RETURN 
 
789
      END 
 
790
C
 
791
 
 
792
 
 
793
 
 
794
      subroutine spin_projection(P,NHEL,IC,TEMP)
 
795
c***********************************************************
 
796
c     This subroutine performs the spin projection.
 
797
c     The relative momentum is recorded in common block
 
798
c     NHEL runs from 1 to nexternal-2, since the sum over 
 
799
c     helicities of the bound state is performed in proj_matrix
 
800
c***********************************************************
 
801
      include 'nexternal.inc'
 
802
      include 'coupl.inc'
 
803
      include 'color_data.inc'
 
804
C
 
805
C     PARAMETERS
 
806
C
 
807
      INTEGER  NEXT_AMP                    ! number of  particles
 
808
      PARAMETER (NEXT_AMP=NEXTERNAL+1)     ! # of external particles before the projection
 
809
C
 
810
C     ARGUMENTS
 
811
C
 
812
      REAL*8 P(0:3,NEXTERNAL)
 
813
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
 
814
c      double precision Ponium(0:3), Pboost(0:3)
 
815
C
 
816
C     LOCAL VARIABLES
 
817
C
 
818
      INTEGER NHEL_AMP(next_amp),IC_AMP(next_amp), L1,L2,J,K
 
819
      DOUBLE COMPLEX JIO(6),TEMP(4,NCOLOR)
 
820
      DOUBLE PRECISION P_AMP(0:3,next_amp),
 
821
     &                  P1(0:3),P2(0:3),PTOT(0:3)   ! MOMENTA OF c AND c~ AND TOTAL MOMENTUM
 
822
      double precision Pthree2,Eq2_quark1,Eq2_quark2
 
823
C
 
824
      COMPLEX*16 JAMP(NCOLOR)
 
825
      DOUBLE PRECISION WEIGHT1, WEIGHT2
 
826
c      DOUBLE COMPLEX PROJAMP(NCOLOR)
 
827
C
 
828
C     global
 
829
C
 
830
      double precision prel(0:3)
 
831
      common /to_pWAVE/prel
 
832
 
 
833
      double precision qmass(2)
 
834
      COMMON/to_qmass/qmass
 
835
 
 
836
c      include 'qmass.inc'
 
837
 
 
838
      Pthree2=P(1,nexternal)**2
 
839
     & +P(2,nexternal)**2+P(3,nexternal)**2
 
840
      if(Pthree2.gt.0.001d0) then
 
841
      write(*,*) 'warning : we are not in the Quarkonium rest frame'
 
842
      endif
 
843
      Eq2_quark1=prel(1)**2+prel(2)**2+prel(3)**2+qmass(1)**2
 
844
      Eq2_quark2=prel(1)**2+prel(2)**2+prel(3)**2+qmass(2)**2
 
845
      PTOT(0)=dsqrt( Eq2_quark1)+dsqrt( Eq2_quark2)
 
846
      PTOT(1)=0d0
 
847
      PTOT(2)=0d0
 
848
      PTOT(3)=0d0
 
849
 
 
850
      weight1=dsqrt(Eq2_quark1)/(dsqrt(Eq2_quark1)+dsqrt(Eq2_quark2))
 
851
      weight2=dsqrt(Eq2_quark2)/(dsqrt(Eq2_quark1)+dsqrt(Eq2_quark2))
 
852
 
 
853
      DO L1=0,3
 
854
      DO L2=1, next_amp-2
 
855
      P_AMP(L1,L2)=P(L1,L2)
 
856
      ENDDO
 
857
C                 in the amplitude
 
858
      P_AMP(L1,next_amp-1)=PTOT(L1)*weight1+prel(L1)
 
859
      P_AMP(L1,next_amp)=PTOT(L1)*weight2-prel(L1)
 
860
      P1(L1)=PTOT(L1)*weight1+prel(L1)
 
861
      P2(L1)=PTOT(L1)*weight2-prel(L1)
 
862
c
 
863
      ENDDO
 
864
C     RECORD THE POLARIZATIONS
 
865
C
 
866
      DO J=1,next_amp-2
 
867
      NHEL_AMP(J)=NHEL(J)
 
868
      IC_AMP(J)=IC(J)
 
869
      ENDDO
 
870
      IC_AMP(next_amp-1)=1
 
871
      IC_AMP(next_amp)=1
 
872
C
 
873
C************************************* 
 
874
C     START THE SUM   
 
875
C************************************* 
 
876
C
 
877
C     INITIALIZATION
 
878
C
 
879
      DO J=1,4,1
 
880
      DO K=1,NCOLOR
 
881
       TEMP(J,K)=(0.0D0,0.0D0)
 
882
      ENDDO
 
883
      ENDDO
 
884
 
 
885
      DO L1=-1,1,2  ! SUM OVER HELICITIES
 
886
      DO L2=-1,1,2
 
887
      NHEL_AMP(next_amp-1)=L1
 
888
      NHEL_AMP(next_amp)=L2  
 
889
      CALL VECTOR0(P1,P2,L1,L2,qmass(1),qmass(2),JIO)! for amp.
 
890
 
 
891
      CALL matrix(P_AMP,NHEL_AMP,IC_AMP,JAMP)         ! for amp.
 
892
C 
 
893
C 
 
894
 
 
895
      DO K=1, NCOLOR    ! INDEX OF JAMP
 
896
          DO J=1,4,1   ! PAY ATTENTION TO THE RANGE : J IN [1,4]
 
897
C                  THE COLOR FACTOR IS INCLUDED IN SQAMP
 
898
          TEMP(J,K)=TEMP(J,K)+JAMP(K)*JIO(J)
 
899
          ENDDO
 
900
      ENDDO
 
901
C 
 
902
C  
 
903
      ENDDO
 
904
      ENDDO
 
905
c
 
906
c     above, J is the spin index, K is the color index
 
907
c
 
908
      return
 
909
      end
 
910
      real*8 function PROJ_H(PP,NHEL,IC)
 
911
C
 
912
C******************************************************************************
 
913
C     THIS FUNCTION PROJECTS THE AMPLITUDE ONTO A P wave CC~ STATE             * 
 
914
C     AND SQUARE IT                                                           *
 
915
C     THE ROUTINE IS UNIVERSAL.                                               *     
 
916
C     MOMENTA MUST BE GIVEN IN THE QUARKONIUM REST FRAME                      *
 
917
C******************************************************************************
 
918
C
 
919
      implicit none
 
920
      include 'nexternal.inc'
 
921
      include 'coupl.inc'
 
922
      include 'color_data.inc'
 
923
C
 
924
C     ARGUMENTS
 
925
C
 
926
      REAL*8 PP(0:3,NEXTERNAL)
 
927
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
 
928
C
 
929
C     LOCAL VARIABLES
 
930
C
 
931
C         1. for the projection
 
932
      DOUBLE PRECISION Pboost(0:3) !,Pboostnd(0:3)
 
933
      REAL*8 P(0:3,NEXTERNAL)
 
934
c      double precision Ponium(0:3),Poniumnd(0:3)
 
935
C
 
936
      DOUBLE COMPLEX TEMP(NCOLOR),
 
937
     & TEMP1(NCOLOR),TEMP2(NCOLOR),TEMP3(NCOLOR)
 
938
C
 
939
      DOUBLE COMPLEX AMP(3,NCOLOR)
 
940
      DOUBLE COMPLEX hel_AMP(NCOLOR),eps(3,-1:1)
 
941
C
 
942
C         2.  to  square the amplitude
 
943
      COMPLEX*16 ZTEMP
 
944
      integer j,k,l,ia,mu
 
945
C
 
946
c--masses
 
947
      DOUBLE PRECISION       qmass(2)
 
948
      COMMON/to_qmass/qmass      
 
949
C
 
950
C     global
 
951
C
 
952
      double precision prel(0:3)
 
953
      common /to_pWAVE/prel
 
954
 
 
955
      double precision pnd(0:3,nexternal),delta
 
956
c      common/modif_p/pnd,delta
 
957
C
 
958
C----
 
959
C  Begin code
 
960
C----
 
961
 
 
962
      delta=0.0001d0      
 
963
c
 
964
c   go to the quarkonium rest frame
 
965
c
 
966
 
 
967
       pboost(0)=pp(0,nexternal)
 
968
c       pboostnd(0)=pnd(0,nexternal)
 
969
       do k=1,3
 
970
       pboost(k)=-pp(k,nexternal)
 
971
c       pboostnd(k)=-pnd(k,nexternal)
 
972
       enddo
 
973
       do j=1,nexternal
 
974
       call boostx(PP(0,j),pboost,P(0,j))
 
975
c       call boostx(Pnd(0,j),pboostnd,Pnd(0,j))
 
976
        do mu=0,3
 
977
        Pnd(mu,j)=P(mu,j)
 
978
        enddo
 
979
        enddo
 
980
 
 
981
        Pnd(0,nexternal)=dsqrt(P(0,nexternal)**2+delta**2)
 
982
 
 
983
c
 
984
c      a) no relative momentum
 
985
c
 
986
       do k=0,3
 
987
       prel(K)=0d0
 
988
       enddo
 
989
c
 
990
 
 
991
       prel(1)=delta/10000d0
 
992
       call spin_singlet_proj(P,NHEL,IC,TEMP)
 
993
 
 
994
c
 
995
c       a) component q_1
 
996
c
 
997
       prel(1)=delta
 
998
 
 
999
       call spin_singlet_proj(Pnd,NHEL,IC,TEMP1)
 
1000
c
 
1001
c       b) component q_2
 
1002
c
 
1003
       prel(1)=0d0
 
1004
       prel(2)=delta
 
1005
 
 
1006
       call spin_singlet_proj(Pnd,NHEL,IC,TEMP2)
 
1007
c
 
1008
c       c) component q_3
 
1009
c
 
1010
       prel(2)=0d0
 
1011
       prel(3)=delta
 
1012
 
 
1013
       call spin_singlet_proj(Pnd,NHEL,IC,TEMP3)
 
1014
c
 
1015
c      computing derivatives
 
1016
c
 
1017
       do k=1,NCOLOR
 
1018
       amp(1,k)=-(temp1(k)-temp(k))/dcmplx(delta)
 
1019
       amp(2,k)=-(temp2(k)-temp(k))/dcmplx(delta)
 
1020
       amp(3,k)=-(temp3(k)-temp(k))/dcmplx(delta)
 
1021
       enddo
 
1022
 
 
1023
C************************************
 
1024
C    Selecting helicity state
 
1025
C    and squaring the amplitude
 
1026
C************************************
 
1027
 
 
1028
      call get_helicity(PP(0,nexternal),eps)
 
1029
 
 
1030
      do L=1,ncolor
 
1031
        hel_amp(l)=(0d0,0d0)
 
1032
        do Ia=1,3
 
1033
          hel_amp(l)=hel_amp(l)+amp(Ia,l)*eps(Ia,nhel(nexternal))
 
1034
        enddo
 
1035
      enddo
 
1036
 
 
1037
      ztemp=(0d0,0d0)
 
1038
      DO K=1,ncolor
 
1039
      DO L=1,ncolor
 
1040
      ztemp=ztemp+hel_amp(K)
 
1041
     &*DCONJG(hel_amp(L)) 
 
1042
     &  *CF(K,L)/DENOM(L)
 
1043
      enddo
 
1044
      enddo
 
1045
 
 
1046
      PROJ_H = DBLE(ztemp)
 
1047
      RETURN 
 
1048
      END 
 
1049
C
 
1050
 
 
1051
 
 
1052
 
 
1053
      subroutine spin_singlet_proj(P,NHEL,IC,TEMP)
 
1054
c***********************************************************
 
1055
c     This subroutine performs the spin projection.
 
1056
c     The relative momentum is recorded in common block
 
1057
c     NHEL runs from 1 to nexternal-1, since the sum over 
 
1058
c     helicities of the bound state is performed in proj_matrix
 
1059
c***********************************************************
 
1060
      include 'nexternal.inc'
 
1061
      include 'coupl.inc'
 
1062
      include 'color_data.inc'
 
1063
C
 
1064
C     PARAMETERS
 
1065
C
 
1066
      INTEGER  NEXT_AMP                    ! number of  particles
 
1067
      PARAMETER (NEXT_AMP=NEXTERNAL+1)     ! # of external particles before the projection
 
1068
C
 
1069
C     ARGUMENTS
 
1070
C
 
1071
      REAL*8 P(0:3,NEXTERNAL)
 
1072
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
 
1073
c      double precision Ponium(0:3), Pboost(0:3)
 
1074
C
 
1075
C     LOCAL VARIABLES
 
1076
C
 
1077
      INTEGER NHEL_AMP(next_amp),IC_AMP(next_amp), L1,L2,J,K
 
1078
      DOUBLE COMPLEX PS,TEMP(NCOLOR)
 
1079
      DOUBLE PRECISION P_AMP(0:3,next_amp),
 
1080
     &                  P1(0:3),P2(0:3),PTOT(0:3)   ! MOMENTA OF c AND c~ AND TOTAL MOMENTUM
 
1081
      double precision Pthree2,Eq2_quark1,Eq2_quark2,weight1,weight2
 
1082
C
 
1083
      COMPLEX*16 JAMP(NCOLOR)
 
1084
c      DOUBLE COMPLEX PROJAMP(NCOLOR)
 
1085
C
 
1086
C     global
 
1087
C
 
1088
 
 
1089
      double precision prel(0:3)
 
1090
      common /to_pWAVE/prel
 
1091
 
 
1092
      DOUBLE PRECISION       qmass(2)
 
1093
      COMMON/to_qmass/qmass
 
1094
 
 
1095
 
 
1096
 
 
1097
      Pthree2=P(1,nexternal)**2+P(2,nexternal)**2
 
1098
     & +P(3,nexternal)**2
 
1099
      if(Pthree2.gt.0.001d0) then
 
1100
      write(*,*) 'warning : we are not in the Quarkonium rest frame'
 
1101
      endif
 
1102
      Eq2_quark1=prel(1)**2+prel(2)**2+prel(3)**2+qmass(1)**2
 
1103
      Eq2_quark2=prel(1)**2+prel(2)**2+prel(3)**2+qmass(2)**2
 
1104
      PTOT(0)=dsqrt(Eq2_quark1)+dsqrt(Eq2_quark2)
 
1105
      PTOT(1)=0d0
 
1106
      PTOT(2)=0d0
 
1107
      PTOT(3)=0d0
 
1108
      DO L1=0,3
 
1109
      DO L2=1, next_amp-2
 
1110
      P_AMP(L1,L2)=P(L1,L2)
 
1111
      ENDDO
 
1112
 
 
1113
      weight1=dsqrt(Eq2_quark1)/(dsqrt(Eq2_quark1)+dsqrt(Eq2_quark2))
 
1114
      weight2=dsqrt(Eq2_quark2)/(dsqrt(Eq2_quark1)+dsqrt(Eq2_quark2))
 
1115
C                 in the amplitude
 
1116
      P_AMP(L1,next_amp-1)=PTOT(L1)*weight1+prel(L1)
 
1117
      P_AMP(L1,next_amp)=PTOT(L1)*weight2-prel(L1)
 
1118
      P1(L1)=PTOT(L1)*weight1+prel(L1)
 
1119
      P2(L1)=PTOT(L1)*weight2-prel(L1)
 
1120
c
 
1121
      ENDDO
 
1122
c
 
1123
c      boost
 
1124
c
 
1125
c      do k=1,next_amp
 
1126
c      call boostx(p_amp(0,k),Ponium,p_amp(0,k))
 
1127
c      enddo
 
1128
      do L1=0,3
 
1129
      P1(L1)=p_amp(L1,next_amp-1)
 
1130
      P2(L1)=p_amp(L1,next_amp)
 
1131
      enddo
 
1132
C
 
1133
C     RECORD THE POLARIZATIONS
 
1134
C
 
1135
      DO J=1,next_amp-2
 
1136
      NHEL_AMP(J)=NHEL(J)
 
1137
      IC_AMP(J)=IC(J)
 
1138
      ENDDO
 
1139
      IC_AMP(next_amp-1)=1
 
1140
      IC_AMP(next_amp)=1
 
1141
C
 
1142
C************************************* 
 
1143
C     START THE SUM   
 
1144
C************************************* 
 
1145
 
 
1146
C
 
1147
C     INITIALIZATION
 
1148
C
 
1149
      DO K=1,NCOLOR
 
1150
       TEMP(K)=(0.0D0,0.0D0)
 
1151
      ENDDO
 
1152
 
 
1153
      DO L1=-1,1,2  ! SUM OVER HELICITIES
 
1154
      DO L2=-1,1,2
 
1155
      NHEL_AMP(next_amp-1)=L1
 
1156
      NHEL_AMP(next_amp)=L2  
 
1157
      CALL pseudoscalar0(P1,P2,L1,L2,
 
1158
     &    qmass(1),qmass(2),PS)              ! for amp.
 
1159
 
 
1160
      CALL matrix(P_AMP,NHEL_AMP,IC_AMP,JAMP)         ! for amp.
 
1161
C 
 
1162
c       
 
1163
C 
 
1164
 
 
1165
      DO K=1, NCOLOR    ! INDEX OF JAMP
 
1166
C                  THE COLOR FACTOR IS INCLUDED IN SQAMP
 
1167
          TEMP(K)=TEMP(K)+JAMP(K)*PS
 
1168
      ENDDO
 
1169
C 
 
1170
      ENDDO
 
1171
      ENDDO
 
1172
c
 
1173
c     above, J is the spin index, K is the color index
 
1174
c
 
1175
      
 
1176
c
 
1177
c     go back to the rest frame
 
1178
c
 
1179
c      do k=1,ncolor
 
1180
c      call boostx(temp(1,k),pboost,temp(1,k))
 
1181
c      enddo
 
1182
c
 
1183
      return
 
1184
      end
 
1185
 
 
1186
 
 
1187
      real*8 function PROJ_PSILEPT(P,NHEL,IC)
 
1188
C
 
1189
C******************************************************************************
 
1190
C     THIS FUNCTION PROJECT THE AMPLITUDE ONTO A CC~(3S1) STATE               * 
 
1191
C     AND SQUARE IT. The quarkonium is decayed into leptons.                  *
 
1192
C     The routine is universal.                                               *
 
1193
C******************************************************************************
 
1194
C
 
1195
      implicit none
 
1196
      include 'nexternal.inc'
 
1197
      include 'coupl.inc'
 
1198
      include 'color_data.inc'
 
1199
C
 
1200
C     PARAMETERS
 
1201
C
 
1202
      INTEGER  NEXT_AMP                    ! initial number of particles
 
1203
      PARAMETER (NEXT_AMP=NEXTERNAL+1)     ! # of external particles before the projection
 
1204
C
 
1205
C     ARGUMENTS
 
1206
C
 
1207
      REAL*8 P(0:3,NEXTERNAL)
 
1208
      INTEGER NHEL(NEXT_amp), IC(NEXTERNAL)    ! number of helicities =nexternal  !
 
1209
C
 
1210
C     LOCAL VARIABLES
 
1211
C
 
1212
C         1. for the projection
 
1213
      INTEGER NHEL_AMP(next_amp),IC_AMP(next_amp), L1,L2,J,K
 
1214
      DOUBLE COMPLEX JIO(6),TEMP(4,NCOLOR), jio_lept(6)
 
1215
      DOUBLE PRECISION P_AMP(0:3,next_amp),
 
1216
     &                  P1(0:3),P2(0:3),PTOT(0:3)   ! MOMENTA OF c AND c~ AND TOTAL MOMENTUM
 
1217
 
 
1218
      DOUBLE PRECISION WEIGHT1, WEIGHT2
 
1219
c      DOUBLE PRECISION EPS(0:3,-1:1)
 
1220
C
 
1221
      COMPLEX*16 JAMP(NCOLOR)
 
1222
      DOUBLE COMPLEX PROJAMP(NCOLOR)
 
1223
C
 
1224
C         3.  to  square the amplitude
 
1225
      COMPLEX*16 ZTEMP
 
1226
c
 
1227
      double precision k1(0:3),k2(0:3),three_k1(3)
 
1228
      common /lepton_stuff/ k1,k2,three_k1
 
1229
C
 
1230
c--masses
 
1231
      DOUBLE PRECISION       qmass(2)
 
1232
      COMMON/to_qmass/qmass      
 
1233
C
 
1234
C     SPECIFYING THE MOMENTA OF C,C~ + THE TOTAL MOMENTUM
 
1235
C
 
1236
 
 
1237
      weight1=qmass(1)/(qmass(1)+qmass(2))
 
1238
      weight2=qmass(2)/(qmass(1)+qmass(2))
 
1239
 
 
1240
      DO L1=0,3
 
1241
      DO L2=1, next_amp-2
 
1242
      P_AMP(L1,L2)=P(L1,L2)
 
1243
      ENDDO
 
1244
      P_AMP(L1,next_amp-1)=P(L1,nexternal)*WEIGHT1
 
1245
      P_AMP(L1,next_amp)=P(L1,nexternal)*WEIGHT2
 
1246
      P1(L1)=P(L1,nexternal)*WEIGHT1
 
1247
      P2(L1)=P(L1,nexternal)*WEIGHT2
 
1248
      PTOT(L1)=P1(L1)+P2(L1)
 
1249
      ENDDO
 
1250
 
 
1251
C
 
1252
C     RECORD THE POLARIZATIONS
 
1253
C
 
1254
      DO J=1,next_amp-2
 
1255
      NHEL_AMP(J)=NHEL(J)
 
1256
      IC_AMP(J)=IC(J)
 
1257
      ENDDO
 
1258
      IC_AMP(next_amp-1)=1
 
1259
      IC_AMP(next_amp)=1
 
1260
C
 
1261
C************************************* 
 
1262
C     START THE SUM   
 
1263
C************************************* 
 
1264
C
 
1265
C     INITIALIZATION
 
1266
C
 
1267
      DO J=1,4,1
 
1268
      DO K=1,NCOLOR
 
1269
       TEMP(J,K)=(0.0D0,0.0D0)
 
1270
      ENDDO
 
1271
      ENDDO
 
1272
 
 
1273
      DO L1=-1,1,2  ! SUM OVER HELICITIES
 
1274
      DO L2=-1,1,2
 
1275
      NHEL_AMP(next_amp-1)=L1
 
1276
      NHEL_AMP(next_amp)=L2  
 
1277
      CALL VECTOR0(P1,P2,L1,L2,qmass(1),qmass(2),JIO)
 
1278
 
 
1279
 
 
1280
      CALL matrix(P_AMP,NHEL_AMP,IC_AMP,JAMP)
 
1281
C 
 
1282
c       
 
1283
      DO K=1, NCOLOR    ! INDEX OF JAMP
 
1284
          DO J=1,4,1   ! PAY ATTENTION TO THE RANGE : J IN [1,4]
 
1285
          TEMP(J,K)=TEMP(J,K)+JAMP(K)*JIO(J)  ! THE COLOR FACTOR IS INCLUDED IN SQAMP
 
1286
          ENDDO
 
1287
      ENDDO
 
1288
C 
 
1289
      ENDDO
 
1290
      ENDDO
 
1291
C       
 
1292
      l1=nhel(next_amp-1)
 
1293
      l2=nhel(next_amp)
 
1294
      call leptonic_current(k1,k2,l1,l2,qmass(1), jio_lept)
 
1295
C
 
1296
      DO K=1,NCOLOR
 
1297
      PROJAMP(K)=TEMP(1,K)*jio_lept(1)-
 
1298
     &  TEMP(2,K)*jio_lept(2)-
 
1299
     & TEMP(3,K)*jio_lept(3)- 
 
1300
     & TEMP(4,K)*jio_lept(4)
 
1301
      ENDDO
 
1302
C     
 
1303
 
 
1304
C************************************ 
 
1305
C    SQUARE THE AMPLITUDE
 
1306
C***********************************
 
1307
C
 
1308
      PROJ_PSILEPT = 0.D0
 
1309
      DO I = 1, NCOLOR
 
1310
          ZTEMP = (0.D0,0.D0)
 
1311
          DO J = 1, NCOLOR
 
1312
              ZTEMP = ZTEMP + CF(J,I)*PROJAMP(J)
 
1313
          ENDDO
 
1314
          PROJ_PSILEPT =PROJ_PSILEPT
 
1315
     &         +DBLE(ZTEMP*DCONJG(PROJAMP(I)))/DENOM(I)
 
1316
      ENDDO
 
1317
c      write(*,*) 'result ', projected_matrix
 
1318
      RETURN 
 
1319
      END 
 
1320
C
 
1321
      subroutine pseudoscalar(P1,P2,LAMBDA1,LAMBDA2,MC,PS)
 
1322
c
 
1323
c This subroutine computes the relevant current for s=0 projection.
 
1324
c
 
1325
c    input : P1 = momentum of c
 
1326
c            P2 = momentum of c~
 
1327
c            LAMBDA1 = helicity of c
 
1328
c            LAMBDA2 = helicity of c~
 
1329
c
 
1330
c    output : PS = the scalar current
 
1331
c           
 
1332
c   
 
1333
      implicit none
 
1334
      double complex fi(6),fo(6),PS, MOMENTUM(2) 
 
1335
      double precision q(0:3),MC,Norm,q2     
 
1336
 
 
1337
      DOUBLE COMPLEX MAT1(1:4,1:4),GAM5(1:4,1:4),MATRIX(1:4,1:4) 
 
1338
c
 
1339
c     "MAT1" IS THE MATRIX   SL(q)+2 E
 
1340
c     "MATRIX" IS THE MATRIX    GAMMA^5 * MAT1
 
1341
c 
 
1342
      DOUBLE PRECISION E
 
1343
      integer J1, J2,K,LAMBDA1,LAMBDA2                         
 
1344
 
 
1345
      DOUBLE PRECISION P1(0:3), P2(0:3)       
 
1346
      double complex cImag, cZero, cOne
 
1347
      parameter( cImag = ( 0.0d0, 1.0d0 ), cZero = ( 0.0d0, 0.0d0 ),
 
1348
     & cOne=(1.0D0,0.0D0) )
 
1349
 
 
1350
      CALL IXXXXX(P1,MC,LAMBDA1,1,FI)     
 
1351
      CALL OXXXXX(P2,MC,LAMBDA2,-1,FO)
 
1352
 
 
1353
      MOMENTUM(1) = -fo(5)+fi(5)
 
1354
      MOMENTUM(2) = -fo(6)+fi(6)
 
1355
 
 
1356
      q(0) = dble( MOMENTUM(1))   
 
1357
      q(1) = dble( MOMENTUM(2))   
 
1358
      q(2) = dimag(MOMENTUM(2))   
 
1359
      q(3) = dimag(MOMENTUM(1))
 
1360
      q2 = q(0)**2-(q(1)**2+q(2)**2+q(3)**2)
 
1361
      E=dsqrt(q2)/2   
 
1362
C
 
1363
C      DEFINITION OF MAT1
 
1364
C 
 
1365
      MAT1(1,1)=2*E
 
1366
      MAT1(1,2)=cZero
 
1367
      MAT1(1,3)=q(0)-q(3)
 
1368
      MAT1(1,4)=-q(1)+cImag*q(2)
 
1369
      MAT1(2,1)=cZero
 
1370
      MAT1(2,2)=2*E
 
1371
      MAT1(2,3)=-q(1)-cImag*q(2)
 
1372
      MAT1(2,4)=q(0)+q(3)
 
1373
      MAT1(3,1)=q(0)+q(3)
 
1374
      MAT1(3,2)=q(1)-cImag*q(2)
 
1375
      MAT1(3,3)=2*E
 
1376
      MAT1(3,4)=cZero
 
1377
      MAT1(4,1)=q(1)+cImag*q(2)
 
1378
      MAT1(4,2)=q(0)-q(3)
 
1379
      MAT1(4,3)=cZero
 
1380
      MAT1(4,4)=2*E 
 
1381
C
 
1382
C     INITIALIZATION OF GAMMA MATRICES
 
1383
C 
 
1384
      DO J1=1,4,1
 
1385
        DO J2=1,4,1
 
1386
          GAM5(J1,J2)=cZero
 
1387
        ENDDO
 
1388
      ENDDO
 
1389
C
 
1390
C     DEFINITION OF GAMMA 5 MATRIX
 
1391
C     CHIRAL REPRESENTATION
 
1392
C:$
 
1393
 
 
1394
c     gamma 5
 
1395
 
 
1396
      GAM5(1,1)=-cOne
 
1397
      GAM5(2,2)=-cOne
 
1398
      GAM5(3,3)=+cOne
 
1399
      GAM5(4,4)=+cOne
 
1400
 
 
1401
c
 
1402
c     DEFINITION OF MAT2 = GAMMA^MU (SL(q)+2E)
 
1403
 
 
1404
      DO J1=1,4,1
 
1405
        DO J2=1,4,1
 
1406
          MATRIX(J1,J2)=cZero
 
1407
          DO K=1,4,1 
 
1408
            MATRIX(J1,J2)=MATRIX(J1,J2)+GAM5(J1,K)*MAT1(K,J2)
 
1409
          ENDDO
 
1410
        ENDDO
 
1411
      ENDDO
 
1412
C
 
1413
 
 
1414
C    NORMALIZATION - does not contain the factor 1/sqrt(E) from the wave function
 
1415
 
 
1416
      Norm=1.0D0/(4*E*DSQRT(2.0D0)*(E+MC))
 
1417
 
 
1418
c    DEFINITION OF THE CURRENT
 
1419
 
 
1420
 
 
1421
      PS=0
 
1422
        DO J1=1,4,1
 
1423
          DO J2 =1,4,1
 
1424
          PS=PS+Norm*FO(J1)*MATRIX(J1,J2)*FI(J2)
 
1425
          ENDDO
 
1426
        ENDDO
 
1427
      return
 
1428
      end
 
1429
      subroutine pseudoscalar0(P1,P2,LAMBDA1,LAMBDA2,M1,M2,PS)
 
1430
c
 
1431
c This subroutine computes the relevant current for s=0 projection.
 
1432
c
 
1433
c    input : P1 = momentum of c
 
1434
c            P2 = momentum of c~
 
1435
c            LAMBDA1 = helicity of c
 
1436
c            LAMBDA2 = helicity of c~
 
1437
c
 
1438
c    output : PS = the scalar current
 
1439
c           
 
1440
c   
 
1441
      implicit none
 
1442
      double complex fi(6),fo(6),PS
 
1443
      double precision M1,M2,N 
 
1444
      integer lambda1,lambda2
 
1445
 
 
1446
      DOUBLE PRECISION P1(0:3), P2(0:3)       
 
1447
      double complex cImag, cZero, cOne
 
1448
      parameter( cImag = ( 0.0d0, 1.0d0 ), cZero = ( 0.0d0, 0.0d0 ),
 
1449
     & cOne=(1.0D0,0.0D0) )
 
1450
 
 
1451
      CALL IXXXXX(P1,M1,LAMBDA1,1,FI)
 
1452
      CALL OXXXXX(P2,M2,LAMBDA2,-1,FO)
 
1453
 
 
1454
      N=dsqrt(8.0d0*M1*M2)
 
1455
 
 
1456
      PS = (-fo(1)*fi(1)-fo(2)*fi(2)+fo(3)*fi(3)+fo(4)*fi(4))/N
 
1457
      return
 
1458
      end
 
1459
      subroutine vector0(P1,P2,LAMBDA1,LAMBDA2,M1,M2,jio)
 
1460
c
 
1461
c This subroutine computes an off-shell vector current from an external
 
1462
c fermion pair.  The vector boson propagator is given in Feynman gauge
 
1463
c for a massless vector and in unitary gauge for a massive vector.
 
1464
c
 
1465
c input:
 
1466
c       complex fi(6)          : flow-in  fermion                   |fi>
 
1467
c       complex fo(6)          : flow-out fermion                   <fo|
 
1468
c       complex gc(2)          : coupling constants                  gvf
 
1469
c       real    vmass          : mass  of OUTPUT vector v
 
1470
c       real    vwidth         : width of OUTPUT vector v
 
1471
c
 
1472
c output:
 
1473
c       complex jio(6)         : vector current          j^mu(<fo|v|fi>)
 
1474
c     
 
1475
      implicit none
 
1476
      double complex fi(6),fo(6),gc(2),jio(6)
 
1477
      double precision q(0:3),q2,p1(0:3),p2(0:3),N
 
1478
      integer lambda1,lambda2
 
1479
      double precision m1,m2
 
1480
      
 
1481
      double complex cImag
 
1482
      parameter( cImag = ( 0.0d0, 1.0d0 ) )
 
1483
 
 
1484
      gc(1)=(1.0d0,0.0d0)
 
1485
      gc(2)=(1.0d0,0.0d0)
 
1486
 
 
1487
      CALL IXXXXX(P1,M1,LAMBDA1,1,FI)     
 
1488
      CALL OXXXXX(P2,M2,LAMBDA2,-1,FO)
 
1489
 
 
1490
      N=dsqrt(8.0d0*M1*M2)
 
1491
      jio(5) = -fo(5)+fi(5)
 
1492
      jio(6) = -fo(6)+fi(6)
 
1493
 
 
1494
      q(0) = dble( jio(5))
 
1495
      q(1) = dble( jio(6))
 
1496
      q(2) = dimag(jio(6))
 
1497
      q(3) = dimag(jio(5))
 
1498
      q2 = q(0)**2-(q(1)**2+q(2)**2+q(3)**2)
 
1499
 
 
1500
 
 
1501
 
 
1502
 
 
1503
            jio(1) = ( gc(1)*( fo(3)*fi(1)+fo(4)*fi(2))
 
1504
     &                +gc(2)*( fo(1)*fi(3)+fo(2)*fi(4)) )/N
 
1505
            jio(2) = (-gc(1)*( fo(3)*fi(2)+fo(4)*fi(1))
 
1506
     &                +gc(2)*( fo(1)*fi(4)+fo(2)*fi(3)) )/N
 
1507
            jio(3) = ( gc(1)*( fo(3)*fi(2)-fo(4)*fi(1))
 
1508
     &                +gc(2)*(-fo(1)*fi(4)+fo(2)*fi(3)))
 
1509
     &               *cImag/N
 
1510
            jio(4) = ( gc(1)*(-fo(3)*fi(1)+fo(4)*fi(2))
 
1511
     &                +gc(2)*( fo(1)*fi(3)-fo(2)*fi(4)) )/N
 
1512
 
 
1513
c 
 
1514
      return
 
1515
      end
 
1516
      subroutine VECTOR(P1,P2,LAMBDA1,LAMBDA2,MQ,jio)
 
1517
c
 
1518
c This subroutine computes the relevant current for s=1 projection.
 
1519
c
 
1520
c    input : P1 = momentum of c
 
1521
c            P2 = momentum of c~
 
1522
c            LAMBDA1 = helicity of c
 
1523
c            LAMBDA2 = helicity of c~
 
1524
c
 
1525
c    output : jio(1->4) = the current,
 
1526
c             jio(5->6) = the total momentum q=p1+p2  
 
1527
c   
 
1528
      implicit none
 
1529
      double complex fi(6),fo(6),jio(6)  
 
1530
      double precision q(0:3),MQ,Norm,q2     
 
1531
 
 
1532
      DOUBLE COMPLEX MAT1(1:4,1:4),GAM(1:4,1:4,0:3),MATRIX(1:4,1:4,0:3) 
 
1533
c
 
1534
c     "MAT1" IS THE MATRIX   SL(q)+2 E
 
1535
c     "MATRIX" IS THE MATRIX    GAMMA^\MU * MAT1
 
1536
c 
 
1537
      DOUBLE PRECISION E
 
1538
      integer J1, J2,K,L,LAMBDA1,LAMBDA2                         
 
1539
      DOUBLE PRECISION P1(0:3), P2(0:3)       
 
1540
      double complex cImag, cZero, cOne
 
1541
      parameter( cImag = ( 0.0d0, 1.0d0 ), cZero = ( 0.0d0, 0.0d0 ),
 
1542
     & cOne=(1.0D0,0.0D0) )
 
1543
 
 
1544
      CALL IXXXXX(P1,MQ,LAMBDA1,1,FI)     
 
1545
      CALL OXXXXX(P2,MQ,LAMBDA2,-1,FO)
 
1546
 
 
1547
      jio(5) = -fo(5)+fi(5)
 
1548
      jio(6) = -fo(6)+fi(6)
 
1549
 
 
1550
      q(0) = dble( jio(5))   
 
1551
      q(1) = dble( jio(6))   
 
1552
      q(2) = dimag(jio(6))   
 
1553
      q(3) = dimag(jio(5))   
 
1554
      q2 = q(0)**2-(q(1)**2+q(2)**2+q(3)**2)
 
1555
      E=dsqrt(q2)/2   
 
1556
C
 
1557
C      DEFINITION OF MAT1
 
1558
C 
 
1559
      MAT1(1,1)=2*E
 
1560
      MAT1(1,2)=cZero
 
1561
      MAT1(1,3)=q(0)-q(3)
 
1562
      MAT1(1,4)=-q(1)+cImag*q(2)
 
1563
      MAT1(2,1)=cZero
 
1564
      MAT1(2,2)=2*E
 
1565
      MAT1(2,3)=-q(1)-cImag*q(2)
 
1566
      MAT1(2,4)=q(0)+q(3)
 
1567
      MAT1(3,1)=q(0)+q(3)
 
1568
      MAT1(3,2)=q(1)-cImag*q(2)
 
1569
      MAT1(3,3)=2*E
 
1570
      MAT1(3,4)=cZero
 
1571
      MAT1(4,1)=q(1)+cImag*q(2)
 
1572
      MAT1(4,2)=q(0)-q(3)
 
1573
      MAT1(4,3)=cZero
 
1574
      MAT1(4,4)=2*E 
 
1575
C
 
1576
C     INITIALIZATION OF GAMMA MATRICES
 
1577
C 
 
1578
      DO J1=1,4,1
 
1579
        DO J2=1,4,1
 
1580
          DO L=0,3,1
 
1581
          GAM(J1,J2,L)=cZero
 
1582
          ENDDO
 
1583
        ENDDO
 
1584
      ENDDO
 
1585
 
 
1586
C
 
1587
C     DEFINITION OF GAMMA MATRICES
 
1588
C     CHIRAL REPRESENTATION
 
1589
C
 
1590
c     gamma 0
 
1591
 
 
1592
      GAM(1,3,0)=cOne
 
1593
      GAM(2,4,0)=cOne
 
1594
      GAM(3,1,0)=cOne
 
1595
      GAM(4,2,0)=cOne
 
1596
 
 
1597
c
 
1598
c     gamma 1
 
1599
c
 
1600
 
 
1601
      GAM(1,4,1)=cOne
 
1602
      GAM(2,3,1)=cOne
 
1603
      GAM(3,2,1)=-cOne
 
1604
      GAM(4,1,1)=-cOne
 
1605
c
 
1606
c     gamma 2
 
1607
c
 
1608
 
 
1609
      GAM(1,4,2)=-cImag
 
1610
      GAM(2,3,2)=cImag
 
1611
      GAM(3,2,2)=cImag
 
1612
      GAM(4,1,2)=-cImag
 
1613
c
 
1614
c
 
1615
c     gamma 3
 
1616
c
 
1617
 
 
1618
      GAM(1,3,3)=cOne
 
1619
      GAM(2,4,3)=-cOne
 
1620
      GAM(3,1,3)=-cOne
 
1621
      GAM(4,2,3)=cOne
 
1622
 
 
1623
c
 
1624
c     DEFINITION OF MAT2 = GAMMA^MU (SL(q)+2E)
 
1625
 
 
1626
      DO L=0,3,1
 
1627
        DO J1=1,4,1
 
1628
          DO J2=1,4,1
 
1629
          MATRIX(J1,J2,L)=cZero
 
1630
            DO K=1,4,1
 
1631
            MATRIX(J1,J2,L)=MATRIX(J1,J2,L)+GAM(J1,K,L)*MAT1(K,J2)
 
1632
            ENDDO
 
1633
          ENDDO
 
1634
        ENDDO
 
1635
      ENDDO
 
1636
C
 
1637
 
 
1638
C    NORMALIZATION - does not contain the factor 1/sqrt(E) from the wave function
 
1639
 
 
1640
      Norm=1.0D0/(4*E*DSQRT(2.0D0)*(E+MQ))
 
1641
 
 
1642
c    DEFINITION OF THE CURRENT
 
1643
 
 
1644
 
 
1645
      DO L=0,3,1
 
1646
        jio(L+1)=0
 
1647
        DO J1=1,4,1
 
1648
          DO J2 =1,4,1
 
1649
          jio(L+1)=jio(L+1)+Norm*FO(J1)*MATRIX(J1,J2,L)*FI(J2)
 
1650
          ENDDO
 
1651
        ENDDO
 
1652
      ENDDO
 
1653
 
 
1654
 
 
1655
      return
 
1656
      end
 
1657
      subroutine vectorBc(P1,P2,LAMBDA1,LAMBDA2,M1,M2,jio)
 
1658
c
 
1659
c This subroutine computes an off-shell vector current from an external
 
1660
c fermion pair.  The vector boson propagator is given in Feynman gauge
 
1661
c for a massless vector and in unitary gauge for a massive vector.
 
1662
c
 
1663
c input:
 
1664
c       complex fi(6)          : flow-in  fermion                   |fi>
 
1665
c       complex fo(6)          : flow-out fermion                   <fo|
 
1666
c       complex gc(2)          : coupling constants                  gvf
 
1667
c       real    vmass          : mass  of OUTPUT vector v
 
1668
c       real    vwidth         : width of OUTPUT vector v
 
1669
c
 
1670
c output:
 
1671
c       complex jio(6)         : vector current          j^mu(<fo|v|fi>)
 
1672
c     
 
1673
      implicit none
 
1674
      double complex fi(6),fo(6),gc(2),jio(6)
 
1675
      double precision q(0:3),q2,p1(0:3),p2(0:3),N
 
1676
      integer lambda1,lambda2
 
1677
      double precision m1,m2
 
1678
      
 
1679
      double precision rZero, rOne
 
1680
      parameter( rZero = 0.0d0, rOne = 1.0d0 )
 
1681
      double complex cImag, cZero
 
1682
      parameter( cImag = ( 0.0d0, 1.0d0 ), cZero = ( 0.0d0, 0.0d0 ) )
 
1683
 
 
1684
      gc(1)=(1.0d0,0.0d0)
 
1685
      gc(2)=(1.0d0,0.0d0)
 
1686
 
 
1687
      CALL IXXXXX(P1,M1,LAMBDA1,1,FI)     
 
1688
      CALL OXXXXX(P2,M2,LAMBDA2,-1,FO)
 
1689
 
 
1690
      N=dsqrt(2d0*2d0*M1*2d0*M2)
 
1691
      jio(5) = -fo(5)+fi(5)
 
1692
      jio(6) = -fo(6)+fi(6)
 
1693
 
 
1694
      q(0) = dble( jio(5))
 
1695
      q(1) = dble( jio(6))
 
1696
      q(2) = dimag(jio(6))
 
1697
      q(3) = dimag(jio(5))
 
1698
      q2 = q(0)**2-(q(1)**2+q(2)**2+q(3)**2)
 
1699
 
 
1700
 
 
1701
 
 
1702
 
 
1703
            jio(1) = ( gc(1)*( fo(3)*fi(1)+fo(4)*fi(2))
 
1704
     &                +gc(2)*( fo(1)*fi(3)+fo(2)*fi(4)) )/N
 
1705
            jio(2) = (-gc(1)*( fo(3)*fi(2)+fo(4)*fi(1))
 
1706
     &                +gc(2)*( fo(1)*fi(4)+fo(2)*fi(3)) )/N
 
1707
            jio(3) = ( gc(1)*( fo(3)*fi(2)-fo(4)*fi(1))
 
1708
     &                +gc(2)*(-fo(1)*fi(4)+fo(2)*fi(3)))
 
1709
     &               *cImag/N
 
1710
            jio(4) = ( gc(1)*(-fo(3)*fi(1)+fo(4)*fi(2))
 
1711
     &                +gc(2)*( fo(1)*fi(3)-fo(2)*fi(4)) )/N
 
1712
 
 
1713
c 
 
1714
      return
 
1715
      end
 
1716
      subroutine leptonic_current(k1,k2,l1,l2,mq, jio)
 
1717
c
 
1718
      implicit none
 
1719
c
 
1720
c     arguments
 
1721
c
 
1722
      double precision k1(0:3),k2(0:3),mq
 
1723
      integer l1,l2
 
1724
      double complex jio(6)
 
1725
c
 
1726
c     local
 
1727
c
 
1728
      double complex fi(6),fo(6),gc(2)
 
1729
      double precision N
 
1730
c
 
1731
      double precision rZero,  pi
 
1732
      parameter( rZero = 0.0d0, PI=3.1415926d0 )
 
1733
      double complex cImag
 
1734
      parameter( cImag = ( 0.0d0, 1.0d0 ))
 
1735
c
 
1736
      N=dsqrt(3.0d0)/(8.0d0*MQ*dsqrt(pi))
 
1737
c
 
1738
      call oxxxxx(k1,rZero,l1,1,FO)
 
1739
      call ixxxxx(k2,rZero,l2,-1,FI)
 
1740
c
 
1741
      gc(1)=(1.0d0,0.0d0)
 
1742
      gc(2)=(1.0d0,0.0d0)
 
1743
c
 
1744
      jio(5) = fo(5)-fi(5)
 
1745
      jio(6) = fo(6)-fi(6)
 
1746
c
 
1747
            jio(1) = ( gc(1)*( fo(3)*fi(1)+fo(4)*fi(2))
 
1748
     &                +gc(2)*( fo(1)*fi(3)+fo(2)*fi(4)) )*N
 
1749
            jio(2) = (-gc(1)*( fo(3)*fi(2)+fo(4)*fi(1))
 
1750
     &                +gc(2)*( fo(1)*fi(4)+fo(2)*fi(3)) )*N
 
1751
            jio(3) = ( gc(1)*( fo(3)*fi(2)-fo(4)*fi(1))
 
1752
     &                +gc(2)*(-fo(1)*fi(4)+fo(2)*fi(3)))
 
1753
     &               *N*cImag
 
1754
            jio(4) = ( gc(1)*(-fo(3)*fi(1)+fo(4)*fi(2))
 
1755
     &                +gc(2)*( fo(1)*fi(3)-fo(2)*fi(4)) )*N
 
1756
c
 
1757
      return
 
1758
      end
 
1759
 
 
1760
      subroutine genmom_lept(P,x,k1,k2,three_k1)
 
1761
C************************************************************************
 
1762
C     GIVEN X(2) AND P SETS UP THE MOMENTUM k1-k2
 
1763
C     AND ALSO SETS UP THE APPROPRIATE
 
1764
C     JACOBIAN VALUE = JAC AND PHASESPACE WEIGTH PSWGT
 
1765
C     X(1) = COSTHETA OF FIRST DECAY psi -> mu mu
 
1766
C     X(2) = PHI OF FIRST DECAY psi -> mu mu
 
1767
C************************************************************************
 
1768
      IMPLICIT NONE
 
1769
C     
 
1770
C     CONSTANTS
 
1771
C
 
1772
      REAL*8     ZERO    ,      PI
 
1773
      PARAMETER (ZERO=0D0,  PI=3.1415926d0 )
 
1774
C
 
1775
C     ARGUMENTS
 
1776
C
 
1777
      REAL*8 P(0:3),k1(0:3),k2(0:3)
 
1778
      REAL*8 three_k1(3)
 
1779
      REAL*8 X(2)
 
1780
C
 
1781
C     LOCAL
 
1782
C
 
1783
      REAL*8 RSH,SH,COSTH,PHI
 
1784
      integer j
 
1785
C
 
1786
C     EXTERNAL
 
1787
C
 
1788
      REAL*8 DOT
 
1789
C
 
1790
C     GLOBAL
 
1791
C
 
1792
      DOUBLE PRECISION              S,X1,X2,PSWGT,JAC
 
1793
      COMMON /PHASESPACE/ S,X1,X2,PSWGT,JAC
 
1794
c
 
1795
      SH = DOT(P,P)
 
1796
      RSH = dsqrt(SH)
 
1797
C
 
1798
C       PICK VALUE OF THETA AND PHI FOR FIRST DECAY SH->M1+S1
 
1799
C
 
1800
        COSTH=-1.D0+2.D0*X(1)
 
1801
        PHI = 2D0*PI*X(2)
 
1802
C
 
1803
C       DETERMINE JACOBIAN FOR THETA AND PHI
 
1804
C     
 
1805
        JAC =  JAC * 4D0*PI
 
1806
C
 
1807
C       CALCULATE COMPONENTS OF MOMENTUM FOUR-VECTORS
 
1808
C       OF THE FINAL STATE MASSLESS PARTONS IN THEIR CM FRAME
 
1809
C
 
1810
        CALL MOM2CX(RSH,zero,zero,COSTH,PHI,k1,k2)
 
1811
c
 
1812
c       record three_k1
 
1813
c
 
1814
        do j=1,3
 
1815
        three_k1(j)=k1(j)
 
1816
        enddo
 
1817
c
 
1818
c       boost the momenta in the lab frame
 
1819
c
 
1820
        call boostx(k1,P,k1)
 
1821
        call boostx(k2,P,k2)
 
1822
 
 
1823
      RETURN
 
1824
      END
 
1825
C
 
1826
C*********************************************************************
 
1827
      SUBROUTINE polariz(P,EPS)
 
1828
C*********************************************************************
 
1829
C 
 
1830
C     ARGUMENTS 
 
1831
      IMPLICIT NONE
 
1832
C 
 
1833
      DOUBLE PRECISION P(0:3),EPS(0:3,-1:1)
 
1834
      
 
1835
C
 
1836
C     LOCAL
 
1837
C 
 
1838
      DOUBLE PRECISION E1(1:3),E2(1:3),E3(1:3)
 
1839
      DOUBLE PRECISION J1(1:3),J2(1:3)
 
1840
      DOUBLE PRECISION ALPHA(3),BETA(2)
 
1841
      INTEGER J
 
1842
      DOUBLE PRECISION SQP, SQNORM_P
 
1843
      DOUBLE PRECISION EPS_SQMIN,EPS_SQMAX,j2_epsmin,norm1,norm2
 
1844
C 
 
1845
C     EXTERNAL
 
1846
C     
 
1847
      DOUBLE PRECISION DOT
 
1848
      EXTERNAL DOT
 
1849
C
 
1850
C      PARAMETER
 
1851
C
 
1852
      DOUBLE PRECISION ZERO
 
1853
      PARAMETER (ZERO=0.0D0)
 
1854
C
 
1855
C     POLARIZATIONS OF J=1 STATE
 
1856
C 
 
1857
      E1(1)=1.0D0
 
1858
      E1(2)=0.0D0
 
1859
      E1(3)=0.0D0
 
1860
C 
 
1861
      E2(1)=0.0D0
 
1862
      E2(2)=1.0D0
 
1863
      E2(3)=0.0D0
 
1864
C 
 
1865
      E3(1)=0.0D0
 
1866
      E3(2)=0.0D0
 
1867
      E3(3)=1.0D0
 
1868
 
 
1869
C
 
1870
C     IF THE STATE IS AT REST, THE SUBROUTINE RETURNS 
 
1871
C     THE THREE VECTORS ABOVE
 
1872
C
 
1873
      IF (P(1).EQ.zero.AND.P(2).eq.zero.AND.P(3).eq.zero) THEN
 
1874
      DO J=1,3,1
 
1875
      EPS(J,-1)=E1(J)
 
1876
      EPS(J,0)=E2(J)
 
1877
      EPS(J,1)=E3(J)
 
1878
      EPS(0,-1)=ZERO
 
1879
      EPS(0,0)=ZERO
 
1880
      EPS(0,1)=ZERO
 
1881
      ENDDO
 
1882
      RETURN
 
1883
      ENDIF
 
1884
C 
 
1885
C     PROJECTIONS
 
1886
C
 
1887
 
 
1888
      ALPHA(1)=P(1)
 
1889
      ALPHA(2)=P(2)
 
1890
      ALPHA(3)=P(3)
 
1891
 
 
1892
      IF (dabs(ALPHA(1)).LT.dabs(ALPHA(3)).AND.dabs(ALPHA(2))
 
1893
     &.LT.dabs(ALPHA(3))) THEN
 
1894
      DO J=1,3,1
 
1895
      J1(J)=E1(J)
 
1896
      J2(J)=E2(J)
 
1897
      ENDDO
 
1898
      BETA(1)=ALPHA(1)
 
1899
      BETA(2)=ALPHA(2)
 
1900
      ELSEIF (dabs(ALPHA(1)).LT.dabs(ALPHA(2)).AND.dabs(ALPHA(3))
 
1901
     & .LT.dabs(ALPHA(2))) THEN
 
1902
      DO J=1,3,1
 
1903
      J1(J)=E1(J)
 
1904
      J2(J)=E3(J)
 
1905
      ENDDO
 
1906
      BETA(1)=ALPHA(1)
 
1907
      BETA(2)=ALPHA(3)
 
1908
      ELSEIF (dabs(ALPHA(2)).LT.dabs(ALPHA(1)).AND.dabs(ALPHA(3)).LT.
 
1909
     & dabs(ALPHA(1))) THEN
 
1910
      DO J=1,3,1
 
1911
      J1(J)=E2(J)
 
1912
      J2(J)=E3(J)
 
1913
      ENDDO 
 
1914
      BETA(1)=ALPHA(2)
 
1915
      BETA(2)=ALPHA(3)
 
1916
      ENDIF
 
1917
C
 
1918
      SQNORM_P=P(1)**2+P(2)**2+P(3)**2
 
1919
C
 
1920
      EPS_SQMIN=0d0
 
1921
      EPS_SQMAX=0d0
 
1922
      j2_epsmin=0d0
 
1923
      DO J=1,3,1
 
1924
      EPS(J,-1)=J1(J)-BETA(1)*P(J)/SQNORM_P
 
1925
C
 
1926
      EPS_SQMIN=EPS_SQMIN+EPS(J,-1)*EPS(J,-1)
 
1927
      j2_epsmin=j2_epsmin+EPS(j,-1)*J2(J) 
 
1928
      ENDDO
 
1929
C
 
1930
      NORM1=DSQRT(EPS_SQMIN)
 
1931
      j2_epsmin=j2_epsmin/norm1
 
1932
C
 
1933
      DO J=1,3
 
1934
      EPS(J,-1)=EPS(J,-1)/norm1
 
1935
      EPS(J,1)=J2(J)-BETA(2)*P(J)/SQNORM_P-j2_epsmin*EPS(J,-1)
 
1936
      EPS_SQMAX=EPS_SQMAX+EPS(J,1)*EPS(J,1)
 
1937
      ENDDO
 
1938
      NORM2=DSQRT(EPS_SQMAX)
 
1939
 
 
1940
      EPS(0,-1)=0d0
 
1941
      EPS(0,1)=0d0
 
1942
 
 
1943
      SQP=P(0)**2-SQNORM_P
 
1944
      DO J=1,3,1
 
1945
      EPS(J,1)=EPS(J,1)/norm2
 
1946
      EPS(J,0)=P(0)*P(J)/DSQRT(SQNORM_P*SQP)
 
1947
      ENDDO 
 
1948
      EPS(0,0)=DSQRT(SQNORM_P/SQP)
 
1949
      RETURN
 
1950
      END
 
1951
 
 
1952
      subroutine get_helicity(P,eps)
 
1953
c
 
1954
c     Given P (momentum in the lab), this subroutine gets the polarization vector 
 
1955
c     Frame: px=py=pz=0 (rest frame of the particle)
 
1956
c     Quantization axis: P
 
1957
 
 
1958
c
 
1959
c     arguments
 
1960
c
 
1961
      double precision P(0:3)
 
1962
      double complex eps(3,-1:1)
 
1963
c
 
1964
c      local:
 
1965
c
 
1966
      double precision pt, pp,sqh
 
1967
c---
 
1968
c Begin code     
 
1969
c---
 
1970
      pt=dsqrt(p(1)**2+p(2)**2)
 
1971
      pp=dsqrt(p(1)**2+p(2)**2+p(3)**2)
 
1972
      sqh=dsqrt(1d0/2d0)
 
1973
c
 
1974
      eps(1,0)=dcmplx(p(1)/pp)
 
1975
      eps(2,0)=dcmplx(p(2)/pp)
 
1976
      eps(3,0)=dcmplx(p(3)/pp)
 
1977
c
 
1978
      eps(1,1)=dcmplx(-p(2)/pt*sqh , p(3)*p(1)/(pt*pp)*sqh)
 
1979
      eps(2,1)=dcmplx(p(1)/pt*sqh , p(3)*p(2)/(pt*pp)*sqh)
 
1980
      eps(3,1)=dcmplx( 0d0 , -pt/pp*sqh )
 
1981
c
 
1982
      eps(1,-1)=dcmplx( p(2)/pt*sqh , p(3)*p(1)/(pt*pp)*sqh)
 
1983
      eps(2,-1)=dcmplx(-p(1)/pt*sqh , p(3)*p(2)/(pt*pp)*sqh)
 
1984
      eps(3,-1)=dcmplx( 0d0 , -pt/pp*sqh )
 
1985
c
 
1986
      return
 
1987
      end
 
1988
 
 
1989
      subroutine get_helicity2(P,eps,m)
 
1990
c
 
1991
c     Given P (momentum in the lab), this subroutine gets the polarization vector
 
1992
c     Frame: lab
 
1993
c     Quantization axis: P
 
1994
 
 
1995
c
 
1996
c     arguments
 
1997
c
 
1998
      double precision P(0:3),m
 
1999
      double complex eps(0:3,-1:1)
 
2000
c
 
2001
c      local:
 
2002
c
 
2003
      double precision pt, pp,sqh
 
2004
c---
 
2005
c Begin code
 
2006
c---
 
2007
      pt=dsqrt(p(1)**2+p(2)**2)
 
2008
      pp=dsqrt(p(1)**2+p(2)**2+p(3)**2)
 
2009
      sqh=dsqrt(1d0/2d0)
 
2010
c 
 
2011
      if (pt.ne.0d0) then
 
2012
      eps(0,0)=dcmplx(pp/m)
 
2013
      eps(1,0)=dcmplx(p(1)/pp*P(0)/m)
 
2014
      eps(2,0)=dcmplx(p(2)/pp*P(0)/m)
 
2015
      eps(3,0)=dcmplx(p(3)/pp*P(0)/m)
 
2016
c
 
2017
      eps(0,1)=dcmplx(0d0)
 
2018
      eps(1,1)=dcmplx(-p(2)/pt*sqh , p(3)*p(1)/(pt*pp)*sqh)
 
2019
      eps(2,1)=dcmplx(p(1)/pt*sqh , p(3)*p(2)/(pt*pp)*sqh)
 
2020
      eps(3,1)=dcmplx( 0d0 , -pt/pp*sqh )
 
2021
c
 
2022
      eps(0,-1)=dcmplx(0d0)
 
2023
      eps(1,-1)=dcmplx( p(2)/pt*sqh , p(3)*p(1)/(pt*pp)*sqh)
 
2024
      eps(2,-1)=dcmplx(-p(1)/pt*sqh , p(3)*p(2)/(pt*pp)*sqh)
 
2025
      eps(3,-1)=dcmplx( 0d0 , -pt/pp*sqh )
 
2026
c
 
2027
      else
 
2028
  
 
2029
      eps(0,0)=dcmplx(p(3)/m)
 
2030
      eps(1,0)=dcmplx(0d0)
 
2031
      eps(2,0)=dcmplx(0d0)
 
2032
      eps(3,0)=dcmplx(P(0)/m)
 
2033
c
 
2034
      eps(0,1)=dcmplx(0d0)
 
2035
      eps(1,1)=dcmplx(-sqh , 0d0)
 
2036
      eps(2,1)=dcmplx( 0d0 , sign(sqh,p(3)))
 
2037
      eps(3,1)=dcmplx( 0d0  )
 
2038
c
 
2039
      eps(0,-1)=dcmplx(0d0)
 
2040
      eps(1,-1)=dcmplx( sqh , 0d0)
 
2041
      eps(2,-1)=dcmplx(0d0 , sign(sqh,p(3)))
 
2042
      eps(3,-1)=dcmplx( 0d0  )
 
2043
      endif
 
2044
 
 
2045
      return
 
2046
      end