1
SUBROUTINE SMATRIX%(proc_id)s(P,ANS)
5
C MadGraph5_aMC@NLO for Madevent Version
7
C Returns amplitude squared summed/avg over colors
9
c for the point in phase space P(0:3,NEXTERNAL)
19
Include 'maxconfigs.inc'
20
Include 'nexternal.inc'
23
PARAMETER ( NCOMB=%(ncomb)d)
25
PARAMETER (NLOOPDIAGRAMS=%(nloopdiags)d)
27
PARAMETER (NLoopFlows=%(nLoopFlows)d)
29
PARAMETER (THEL=NCOMB)
30
INTEGER Hel_Average_Factor
31
PARAMETER (Hel_Average_Factor=%(hel_avg_factor)d)
35
REAL*8 P(0:3,NEXTERNAL),ANS
37
c global (due to reading writting)
39
LOGICAL GOODHEL(NCOMB)
41
common/BLOCK_GOODHEL/NTRY,GOODHEL
45
INTEGER NHEL(NEXTERNAL,NCOMB),NTRY
46
REAL*8 T,MATRIX%(proc_id)s
47
REAL*8 R,SUMHEL,TS(NCOMB)
49
INTEGER IPROC,JC(NEXTERNAL),II
50
LOGICAL GOODHEL(NCOMB)
51
REAL*8 HWGT, XTOT, XTRY, XREJ, XR, YFRAC(0:NCOMB)
52
INTEGER IDUM, NGOOD, IGOOD(NCOMB), JHEL, J, JJ
56
REAL*8 , ALLOCATABLE :: ANS_ML(:,:)
57
REAL*8 , ALLOCATABLE :: PREC_FOUND(:)
59
INTEGER ML_ANS_DIMENSION
60
C Prec_found ans ans_ml are saved so that they don't have to be re-allocated everytime
61
SAVE ML_ANS_DIMENSION,ANS_ML,PREC_FOUND
64
DATA XTRY, XREJ, NGOOD /0,0,0/
67
DATA INIT_MADLOOP/.FALSE./
71
LOGICAL FORCE_ML_HELICITY_SUM
72
common/to_ML_control/FORCE_ML_HELICITY_SUM
74
DOUBLE PRECISION AMP2(MAXAMPS), JAMP2(0:MAXFLOW)
75
COMMON/TO_AMPS/ AMP2, JAMP2
77
CHARACTER*101 HEL_BUFF
78
COMMON/TO_HELICITY/ HEL_BUFF
80
%(real_dp_format)s ML_JAMP2(nLoopFlows)
81
common/%(ml_prefix)sJAMP2/ML_JAMP2
83
%(real_dp_format)s ML_AMP2(NLOOPDIAGRAMS)
84
common/%(ml_prefix)sAMP2/ML_AMP2
87
COMMON/TO_POLARIZATION/ POL
91
COMMON/TO_MATRIX/ISUM_HEL, MULTI_CHANNEL
92
%(define_iconfigs_lines)s
93
DATA XTRY, XREJ, NGOOD /0,0,0/
94
SAVE YFRAC, IGOOD, JHEL
102
IF(NLOOPDIAGRAMS.GT.MAXAMPS) THEN
103
STOP 'MAXAMPS smaller than to NLOOPDIAGRAMS in matrix.f'
106
IF(.NOT.INIT_MADLOOP) THEN
107
INIT_MADLOOP = .TRUE.
108
C So that TIR integrals can be reused accross helicity configuration.
109
C Notice that ITR integrals for rotated PS points (i.e. which is part of MadLoop's default stabiliy test procedure), will not work. This is why one should always keep the MadLoop parameter 'NRotations_DP' set to 0.
110
CALL %(ml_prefix)sSET_AUTOMATIC_TIR_CACHE_CLEARING(.FALSE.)
111
CALL %(ml_prefix)sGET_ANSWER_DIMENSION(ML_ANS_DIMENSION)
112
ALLOCATE(ANS_ML(0:3,0:ML_ANS_DIMENSION))
113
ALLOCATE(PREC_FOUND(0:ML_ANS_DIMENSION))
115
C But then we must clear the cache by hand at the beginning of the computation of each new PS point
116
CALL %(ml_prefix)sCLEAR_TIR_CACHE()
123
IF (multi_channel) THEN
133
WRITE(HEL_BUFF,'(20I5)') (0,I=1,NEXTERNAL)
139
C Notice that when forcing helicity sum directly in MadLoop then one doesn't have access to the weights of individual Hel. Configs so that they cannot be specified in the event file. We turned this option off by default.
140
FORCE_ML_Helicity_Sum = (ISUM_HEL.eq.0.AND..False.)
141
if (FORCE_ML_Helicity_Sum)then
142
C Of course this option can also not be used for polarized beams.
144
if (POL(JJ).NE.1d0) FORCE_ML_Helicity_Sum = .false.
148
IF(Force_ML_Helicity_Sum )THEN
150
CALL %(ml_prefix)sSLOOPMATRIX_THRES(P,ANS_ML,-1.0d0,PREC_FOUND,ML_RET_CODE)
154
IF (ISUM_HEL.EQ.0.or.(DS_get_dim_status('Helicity').eq.0)) THEN
156
IF (GOODHEL(I) .OR. NTRY .LE. MAXTRIES.OR.(ISUM_HEL.NE.0)) THEN
158
CALL %(ml_prefix)sSLOOPMATRIXHEL_THRES(P,I,ANS_ML,-1.0d0,PREC_FOUND,ML_RET_CODE)
161
JAMP2(JJ) = JAMP2(JJ) + ML_JAMP2(JJ)
163
DO JJ=1,NLOOPDIAGRAMS
164
AMP2(JJ) = AMP2(JJ)+ ML_AMP2(JJ)
168
IF(POL(JJ).NE.1d0.AND.NHEL(JJ,I).EQ.INT(SIGN(1d0,POL(JJ)))) THEN
170
ELSE IF(POL(JJ).NE.1d0)THEN
171
T=T*(2d0-ABS(POL(JJ)))
174
IF (ISUM_HEL.NE.0) then
175
call DS_add_entry('Helicity',I,T)
181
IF (ISUM_HEL.NE.0) then
182
! We set HEL_PICKED to -1 here so that later on, the call to DS_add_point in dsample.f does not add anything to the grid since it was already done here.
184
! For safety, hardset the helicity sampling jacobian to 0.0d0 to make sure it is not .
186
IF(DS_get_dim_status('Helicity').eq.1) then
187
! If we finished the initialization we can update the grid so as to start sampling over it.
188
! However the grid will now be filled by dsample with different kind of weights (including pdf, flux, etc...) so by setting the grid_mode of the reference grid to 'initialization' we make sure it will be overwritten (as opposed to 'combined') by the running grid at the next update.
189
CALL DS_UPDATE_GRID('Helicity')
190
CALL DS_SET_GRID_MODE('Helicity','init')
191
call reset_cumulative_variable() ! avoid biais of the initialization
195
IF(NTRY.LE.MAXTRIES)THEN
197
IF (.NOT.GOODHEL(I) .AND. (TS(I).GT.ANS*LIMHEL/NCOMB)) THEN
201
print *,'Adding good helicity ',I,TS(I)/ANS
205
IF(NTRY.EQ.MAXTRIES)THEN
206
ISUM_HEL=MIN(ISUM_HEL,NGOOD)
208
IF(NTRY(IMIRROR).EQ.(MAXTRIES+1)) THEN
209
call reset_cumulative_variable() ! avoid biais of the initialization
213
ELSE !RANDOM HELICITY
214
C The helicity configuration was chosen already by genps and put in a common block defined in genps.inc.
217
CALL %(ml_prefix)sSLOOPMATRIXHEL_THRES(P,I,ANS_ML,-1.0d0,PREC_FOUND,ML_RET_CODE)
220
JAMP2(JJ) = JAMP2(JJ) + ML_JAMP2(JJ)
222
DO JJ=1,NLOOPDIAGRAMS
223
AMP2(JJ) = AMP2(JJ)+ ML_AMP2(JJ)
227
IF(POL(JJ).NE.1d0.AND.NHEL(JJ,I).EQ.INT(SIGN(1d0,POL(JJ)))) THEN
229
ELSE IF(POL(JJ).NE.1d0)THEN
230
T=T*(2d0-ABS(POL(JJ)))
233
c Always one helicity at a time
235
c Include the Jacobian from helicity sampling
236
ANS = ANS * hel_jacobian
237
WRITE(HEL_BUFF,'(20i5)')(NHEL(II,I),II=1,NEXTERNAL)
239
IF (ISUM_HEL .NE. 1.or.(HEL_PICKED.eq.-1)) THEN
243
SUMHEL=SUMHEL+DABS(TS(I))
245
WRITE(HEL_BUFF,'(20i5)')(NHEL(II,I),II=1,NEXTERNAL)
254
IF (MULTI_CHANNEL) THEN
259
IF (XTOT.NE.0D0) THEN
265
IF(.not.Force_ML_Helicity_Sum )THEN
266
ANS = ANS/ Hel_Average_Factor
269
C Amplitude(s) for diagram number %(n_tot_diags)d
270
C This last line is a tag do not remove it.
added
removed
madgraph/iolibs/template_files/matrix_loop_induced_madevent.inc
