← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.x_eva_v2

« back to all changes in this revision

Viewing changes to Template/LO/SubProcesses/Dummy_ftc_src/dummy_fct_EVA_3A_matching.f

  • Committer: Richard Ruiz
  • Date: 2021-10-15 05:46:45 UTC
  • Revision ID: richard.physics@gmail.com-20211015054645-gbhs3h8952961d65
Adding dummy_ftc repository to Template/LO/SubProcesses

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/LO/SubProcesses/Dummy_ftc_src/dummy_fct_EVA_3A_matching.f
 
1
C**************************************************************************
 
2
c     dummy_fct_EVA_3A_matching.f
 
3
c     2021 October (R. Ruiz)
 
4
c     Imposes phase space cuts on AAA final-state and on (mu-v_mu) momentum 
 
5
c     transfer (q) or transverse momentum of v_mu (pT) for matrix element 
 
6
c     matching of the processes:
 
7
c     1. e mu > v_e v_mu AAA 
 
8
c     2. e W > v_e AAA
 
9
c     pT/q cut are set to factorization scale (dsqrt_q2fact2) in run_card.dat. 
 
10
c     For further details, see Constantini et al [arXiv:2110.]
 
11
C**************************************************************************
 
12
 
 
13
      logical FUNCTION dummy_cuts(P)
 
14
C**************************************************************************
 
15
C     INPUT:
 
16
C     P(0:3,1)           MOMENTUM OF INCOMING PARTON
 
17
C     P(0:3,2)           MOMENTUM OF INCOMING PARTON
 
18
C     P(0:3,3)           MOMENTUM OF ...
 
19
C     ALL MOMENTA ARE IN THE REST FRAME!!
 
20
C     COMMON/JETCUTS/   CUTS ON JETS
 
21
C     OUTPUT:
 
22
C     TRUE IF EVENTS PASSES ALL CUTS LISTED
 
23
C**************************************************************************
 
24
      IMPLICIT NONE
 
25
c     
 
26
c     Constants
 
27
c     
 
28
      include 'genps.inc'
 
29
      include 'nexternal.inc'
 
30
      include 'run.inc'
 
31
      include 'maxamps.inc'
 
32
      integer i,j,k
 
33
      integer idup(nexternal,maxproc,maxsproc)
 
34
      integer mothup(2,nexternal)
 
35
      integer icolup(2,nexternal,maxflow,maxsproc)
 
36
      include 'leshouche.inc'
 
37
C     
 
38
C     ARGUMENTS
 
39
C     
 
40
      REAL*8 P(0:3,nexternal)
 
41
C     
 
42
C     PARAMETERS
 
43
C     
 
44
      real*8 PI
 
45
      parameter( PI = 3.14159265358979323846d0 )
 
46
 
 
47
      integer ff,ii,kk,nrPho
 
48
      double precision evoThres2,evoScale2 ! neutrino
 
49
      double precision pAAA(0:3),minM2,minPtX,maxEta ! photon
 
50
      double precision rap,pt,r2,dot,SumDot
 
51
      external rap,pt,r2,dot,SumDot
 
52
 
 
53
      integer tmpPID1,tmpPID2
 
54
      integer ipdg(nexternal)
 
55
c     
 
56
c     particle identification
 
57
c     
 
58
      logical  doEVACuts,doEVAMatch
 
59
      LOGICAL  IS_A_J(NEXTERNAL),IS_A_L(NEXTERNAL)
 
60
      LOGICAL  IS_A_B(NEXTERNAL),IS_A_A(NEXTERNAL),IS_A_ONIUM(NEXTERNAL)
 
61
      LOGICAL  IS_A_NU(NEXTERNAL),IS_HEAVY(NEXTERNAL)
 
62
      logical  do_cuts(nexternal)
 
63
      COMMON /TO_SPECISA/IS_A_J,IS_A_A,IS_A_L,IS_A_B,IS_A_NU,IS_HEAVY,
 
64
     .     IS_A_ONIUM, do_cuts
 
65
 
 
66
 
 
67
c     (re)set thresholds
 
68
      doEVACuts=.true.
 
69
      doEVAMatch=.true.         ! set to .true.(.false.) for emu (eW) scattering
 
70
      nrPho=0
 
71
      maxEta=3.d0      
 
72
      minPtX=80.d0
 
73
      minM2=(1d3)**2
 
74
      evoScale2 = 0d0
 
75
      evoThres2 = q2fact(2)
 
76
 
 
77
 
 
78
c     for each external:
 
79
      do ff=nincoming+1,nexternal
 
80
 
 
81
c     Cut on evolution variable:
 
82
c     1. Find v_mu (exists only if full ME)
 
83
c     2. Build evolution variable (0=q^2 or 1=pT^2)
 
84
c     3. Cut on evolution variable
 
85
         tmpPID1 = abs(idup(ff,1,1))
 
86
         if(tmpPID1.eq.14.and.doEVAMatch) then
 
87
            if(ievo_eva.eq.1) then
 
88
               evoScale2 = pT(p(0,ff))**2
 
89
            else
 
90
               do ii=1,nincoming
 
91
                  tmpPID2=abs(idup(ii,1,1))
 
92
                  if(tmpPID2.eq.13) then
 
93
                     evoScale2 = abs(SumDot(p(0,ff),p(0,ii),-1d0))
 
94
                     exit
 
95
                  endif
 
96
               enddo
 
97
            endif
 
98
            if(evoScale2.lt.evoThres2) then
 
99
               dummy_cuts=.false.
 
100
               return
 
101
            endif
 
102
         endif
 
103
 
 
104
 
 
105
c     Cut on photons:
 
106
c     1. Find photons and add momentum
 
107
c     2. Once three photons found: impose cut on invariant mass
 
108
c     3. pT/eta cuts are redundant/sanity check
 
109
         if(IS_A_A(ff).and.doEVACuts) then
 
110
            nrPho=nrPho+1
 
111
            if(nrPho.lt.2) then ! firstTime                       
 
112
               do kk=0,3
 
113
                  pAAA(kk) = p(kk,ff) ! p(4-momentum,list-ID)
 
114
               enddo
 
115
            else
 
116
               do kk=0,3
 
117
                  pAAA(kk) = pAAA(kk) + p(kk,ff) ! p(4-momentum,list-ID)
 
118
               enddo
 
119
            endif
 
120
c     check M(AAA)            
 
121
            if(nrPho.gt.2) then
 
122
               if(dot(pAAA,pAAA).lt.minM2) then
 
123
                  dummy_cuts=.false.
 
124
                  return
 
125
               endif
 
126
            endif
 
127
c     check eta of A        
 
128
            if(abs(rap(p(0,ff))).gt.maxEta) then
 
129
               dummy_cuts=.false.
 
130
               return
 
131
            endif
 
132
c     check pT of A
 
133
            if(pt(p(0,ff)).lt.minPtX) then
 
134
               dummy_cuts=.false.
 
135
               return
 
136
            endif            
 
137
         endif         
 
138
      enddo
 
139
 
 
140
      dummy_cuts=.true.
 
141
      return
 
142
      end
 
143
 
 
144
      subroutine get_dummy_x1(sjac, X1, R, pbeam1, pbeam2, stot, shat)
 
145
      implicit none
 
146
      include 'maxparticles.inc'
 
147
      include 'run.inc'
 
148
c     include 'genps.inc'
 
149
      double precision sjac     ! jacobian. should be updated not reinit
 
150
      double precision X1       ! bjorken X. output
 
151
      double precision R        ! random value after grid transfrormation. between 0 and 1
 
152
      double precision pbeam1(0:3) ! momentum of the first beam (input and/or output)
 
153
      double precision pbeam2(0:3) ! momentum of the second beam (input and/or output)
 
154
      double precision stot     ! total energy  (input and /or output)
 
155
      double precision shat     ! output
 
156
 
 
157
c     global variable to set (or not)
 
158
      double precision cm_rap
 
159
      logical set_cm_rap
 
160
      common/to_cm_rap/set_cm_rap,cm_rap
 
161
      
 
162
      set_cm_rap=.false.        ! then cm_rap will be set as .5d0*dlog(xbk(1)*ebeam(1)/(xbk(2)*ebeam(2)))
 
163
!     ebeam(1) and ebeam(2) are defined here thanks to 'run.inc'
 
164
      shat = x1*ebeam(1)*ebeam(2)
 
165
      return 
 
166
      end
 
167
 
 
168
      subroutine get_dummy_x1_x2(sjac, X, R, pbeam1, pbeam2, stot,shat)
 
169
      implicit none
 
170
      include 'maxparticles.inc'
 
171
      include 'run.inc'
 
172
c     include 'genps.inc'
 
173
      double precision sjac     ! jacobian. should be updated not reinit
 
174
      double precision X(2)     ! bjorken X. output
 
175
      double precision R(2)     ! random value after grid transfrormation. between 0 and 1
 
176
      double precision pbeam1(0:3) ! momentum of the first beam
 
177
      double precision pbeam2(0:3) ! momentum of the second beam
 
178
      double precision stot     ! total energy
 
179
      double precision shat     ! output
 
180
 
 
181
c     global variable to set (or not)
 
182
      double precision cm_rap
 
183
      logical set_cm_rap
 
184
      common/to_cm_rap/set_cm_rap,cm_rap
 
185
      
 
186
      set_cm_rap=.false.        ! then cm_rap will be set as .5d0*dlog(xbk(1)*ebeam(1)/(xbk(2)*ebeam(2)))
 
187
!     ebeam(1) and ebeam(2) are defined here thanks to 'run.inc'
 
188
      shat = x(1)*x(2)*ebeam(1)*ebeam(2)
 
189
      return 
 
190
      end
 
191
 
 
192
 
 
193
      logical  function dummy_boostframe()
 
194
      implicit none
 
195
c     
 
196
c     
 
197
      dummy_boostframe = .false.
 
198
      return
 
199
      end
 
200