Merge lp:~maddevelopers/mg5amcnlo/3.0.2-alpha0 into lp:~maddevelopers/mg5amcnlo/3.2.1
- 3.0.2-alpha0
- Merge into 3.2.1
Status: | Merged |
---|---|
Merged at revision: | 973 |
Proposed branch: | lp:~maddevelopers/mg5amcnlo/3.0.2-alpha0 |
Merge into: | lp:~maddevelopers/mg5amcnlo/3.2.1 |
Diff against target: |
2463 lines (+1387/-325) 21 files modified
Template/NLO/Cards/run_card.dat (+4/-1) Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f (+290/-19) Template/NLO/SubProcesses/BinothLHA.f (+37/-3) Template/NLO/SubProcesses/chooser_functions.f (+9/-0) Template/NLO/SubProcesses/cuts.f (+537/-277) Template/NLO/SubProcesses/fks_singular.f (+31/-7) Template/NLO/SubProcesses/makefile_fks_dir (+1/-0) Template/NLO/SubProcesses/recmom.f (+64/-8) UpdateNotes.txt (+2/-1) madgraph/core/diagram_generation.py (+19/-3) madgraph/fks/fks_base.py (+8/-0) madgraph/fks/fks_common.py (+12/-3) madgraph/fks/fks_helas_objects.py (+2/-0) madgraph/fks/fks_tag.py (+75/-0) madgraph/interface/madgraph_interface.py (+37/-3) madgraph/iolibs/export_fks.py (+86/-0) madgraph/iolibs/template_files/fks_info.inc (+7/-0) madgraph/iolibs/template_files/rescale_alpha_tagged.inc (+40/-0) madgraph/various/banner.py (+11/-0) tests/acceptance_tests/test_cmd_amcatnlo.py (+35/-0) tests/unit_tests/fks/test_fks_taggedphotons.py (+80/-0) |
To merge this branch: | bzr merge lp:~maddevelopers/mg5amcnlo/3.0.2-alpha0 |
Related bugs: |
Reviewer | Review Type | Date Requested | Status |
---|---|---|---|
Rikkert Frederix | Approve | ||
Olivier Mattelaer | Approve | ||
Review via email: mp+408723@code.launchpad.net |
Commit message
Description of the change
This merge includes the possibility to use tagged particles in the final state, which is particularly relevant for EW corrections with tagged photons (see 2106.02059).
Note that cuts.f has been deeply refactored.
- 1000. By marco zaro
-
removed process folders
Rikkert Frederix (frederix) wrote : | # |
- 1001. By marco zaro
-
merged with 3.2.1
marco zaro (marco-zaro) wrote : | # |
Hi Rikkert,
thanks.
> On 16 Sep 2021, at 16:22, Rikkert Frederix <email address hidden> wrote:
>
> Hi Marco,
>
> I'll have a closer look next week. But for now:
>
> - There is a text conflict in the UpdateNotes.txt
I solved it, putting the changes for 3.1.2 together with 3.2.1, since the former was never released
> - Why do you add the new models? Shouldn't they part of FeynRules webpage so that they are simply downloaded when needed?
I will discuss with Olivier and put the models on the web repository
Cheers,
Marco
>
> Best,
> Rikkert
>
> --
> https:/
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.2.1.
>
Olivier Mattelaer (olivier-mattelaer) wrote : | # |
Hi,
I have given a quick look and this sounds very good.
- I agree with Rikkert that it is better to put the model online (maybe not on the Feynrules website but on our madgraph website) the branch lp:~maddevelopers/mg5amcnlo/UFOmodel_db is designed for that actually.
- I'm worried about the re-introduction of the fortran function SORTTI and SORTTC.
If I remember correctly, such function had some issue with GCC10, did you check.
I do not see in the diff any kind of test (maybe it is just truncated). Can you add at least on parralel test such that I can run it on various GCC version and check that they are no issue on that.
Did you add unittest for the new syntax? if not I will likely add those when running some tests.
Cheers,
Olivier
- 1002. By marco zaro
-
fix in cuts.f to remove sortti/sorttc; unit+acceptance test added; models removed and added to the UFOmodel branch
marco zaro (marco-zaro) wrote : | # |
Hello both
thanks for the reviews.
I have removed the models rom the branch and added to the UFOmodel one. I should have taken care of all the other comments in revision 1002.
Cheers,
Marco
> On 17 Sep 2021, at 09:54, Olivier Mattelaer <email address hidden> wrote:
>
> Hi,
>
> I have given a quick look and this sounds very good.
>
> - I agree with Rikkert that it is better to put the model online (maybe not on the Feynrules website but on our madgraph website) the branch lp:~maddevelopers/mg5amcnlo/UFOmodel_db is designed for that actually.
>
> - I'm worried about the re-introduction of the fortran function SORTTI and SORTTC.
> If I remember correctly, such function had some issue with GCC10, did you check.
>
> I do not see in the diff any kind of test (maybe it is just truncated). Can you add at least on parralel test such that I can run it on various GCC version and check that they are no issue on that.
>
> Did you add unittest for the new syntax? if not I will likely add those when running some tests.
>
> Cheers,
>
> Olivier
> --
> https:/
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.2.1.
>
Rikkert Frederix (frederix) wrote : | # |
Hi Marco,
I've done the following test with this branch:
import model loop_qcd_
generate p p > !a! !a! [real=QCD] ; output ; launch -p -f
and it crashes stating that there are "no integration channels". Looking in the relevant gensym.log, there is "ERROR, too many photons 2". I thought that the above process should have worked, no? Have I misunderstood something here?
Best,
Rikkert
- 1003. By marco zaro
-
cuts.f updated so that recombine_momenta won't crash for too many photons (thanks Ioannis Tsinikos); the gmu<>alpha scheme conversion is written dinamically, to ensure compatibility with processes without tagged photons
marco zaro (marco-zaro) wrote : | # |
Hi Rikkert,
Indeed, Ioannis had produces a modified version of cuts.f that makes disappear your error. I have updated cuts.f in the branch accordingly, and fixed some other things in the code.
Cheers,
Marco
> On 21 Sep 2021, at 09:56, Rikkert Frederix <email address hidden> wrote:
>
> Hi Marco,
>
> I've done the following test with this branch:
>
> import model loop_qcd_
> generate p p > !a! !a! [real=QCD] ; output ; launch -p -f
>
> and it crashes stating that there are "no integration channels". Looking in the relevant gensym.log, there is "ERROR, too many photons 2". I thought that the above process should have worked, no? Have I misunderstood something here?
>
> Best,
> Rikkert
>
>
> --
> https:/
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.2.1z
Olivier Mattelaer (olivier-mattelaer) wrote : | # |
This is very good and basically ready to be merged.
Some minor point.
1) !z! works ... should not be this forbidden?
2) can we use tagged photon with QCD correction only?
-- In the commit I did, I have added unittest that crash related to the above command, you can remove those test if those are fine--
3) 2!a! was making the code to crash but !2a! was working fine. I have updated the parser to support both.
4) I have spotted a typo in banner.py (see the comment within the diff for that)
5) I have loose track of all the change in cuts.f but did not spot anything obviously wrong. If you are not confident enough about that part maybe I should go deeper and track that we did not "loose" any previous cut but sounds like it.
Thanks a lot and sorry for being slow on this review
- 1004. By olivier-mattelaer
-
add unittest --not all passing-- and update the parser to allow duplicate tagged particcle
- 1005. By marco zaro
-
fix in banner.py according to Olivier's comment
marco zaro (marco-zaro) wrote : | # |
Hi Olivier,
thanks for your comments
> On 14 Oct 2021, at 16:26, Olivier Mattelaer <email address hidden> wrote:
>
> Review: Approve
>
> This is very good and basically ready to be merged.
>
> Some minor point.
> 1) !z! works ... should not be this forbidden?
yes, but it should do nothing
> 2) can we use tagged photon with QCD correction only?
> -- In the commit I did, I have added unittest that crash related to the above command, you can remove those test if those are fine—
In principle the syntax should work, giving identical results with the case without tagging
> 3) 2!a! was making the code to crash but !2a! was working fine. I have updated the parser to support both.
I was not aware of the fact that one can use the numbers to specify multiple times the same particle
> 4) I have spotted a typo in banner.py (see the comment within the diff for that)
I have fixed it, please have a look
> 5) I have loose track of all the change in cuts.f but did not spot anything obviously wrong. If you are not confident enough about that part maybe I should go deeper and track that we did not "loose" any previous cut but sounds like it.
>
Let us see what Rik says here…
Cheers,
Marco
> Thanks a lot and sorry for being slow on this review
>
> Diff comments:
>
>>
>> === modified file 'madgraph/
>> --- madgraph/
>> +++ madgraph/
>> @@ -4552,11 +4553,20 @@
>> if proc_characteri
>> #remove all cut
>> self.remove_
>> +
>> + # check for tagged photons
>> + tagged_particles = set()
>>
>> # Check if need matching
>> min_particle = 99
>> max_particle = 0
>> for proc in proc_def:
>> + for leg in proc['legs']:
>> + if leg['is_tagged']:
>> + tagged_
>> +
>> + if 22 in tagged_particles:
>> + self['gamma_is_j'] = False
>> min_particle = min(len(
>> max_particle = max(len(
>
> The min_particle/
> for proc in proc_def: loop
>
>> matching = False
>
>
> --
> https:/
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.2.1.
>
Rikkert Frederix (frederix) wrote : | # |
Hi Marco,
I think this is ready for merging.
best,
Rikkert
Preview Diff
1 | === modified file 'Template/NLO/Cards/run_card.dat' |
2 | --- Template/NLO/Cards/run_card.dat 2021-04-05 16:07:15 +0000 |
3 | +++ Template/NLO/Cards/run_card.dat 2021-10-15 11:42:16 +0000 |
4 | @@ -129,7 +129,10 @@ |
5 | %(bwcutoff)s = bwcutoff |
6 | #*********************************************************************** |
7 | # Cuts on the jets. Jet clustering is performed by FastJet. * |
8 | -# - If gamma_is_j, photons are also clustered * |
9 | +# - If gamma_is_j, photons are also clustered with jets. * |
10 | +# Otherwise, they will be treated as tagged particles and photon * |
11 | +# isolation will be applied. Note that photons in the real emission * |
12 | +# will always be clustered with QCD partons. * |
13 | # - When matching to a parton shower, these generation cuts should be * |
14 | # considerably softer than the analysis cuts. * |
15 | # - More specific cuts can be specified in SubProcesses/cuts.f * |
16 | |
17 | === modified file 'Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f' |
18 | --- Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f 2018-10-26 09:55:25 +0000 |
19 | +++ Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f 2021-10-15 11:42:16 +0000 |
20 | @@ -5,8 +5,12 @@ |
21 | integer nwgt |
22 | character*(*) weights_info(*) |
23 | integer i,l |
24 | - character*6 cc(2) |
25 | - data cc/'|T@NLO','|T@LO '/ |
26 | +c character*6 cc(2) |
27 | +c data cc/'|T@NLO','|T@LO '/ |
28 | + character*13 cc(9) |
29 | + data cc/ ' |T@NLO ',' |T@LO ',' |T@LO1 ' |
30 | + $ ,' |T@LO2 ',' |T@LO3 ',' |T@NLO1 ' |
31 | + $ ,' |T@NLO2 ',' |T@NLO3 ',' |T@NLO4 '/ |
32 | |
33 | c Also specific perturbative orders can be directly plotted, adding for examples |
34 | c the following further entries in the variable data |
35 | @@ -18,8 +22,8 @@ |
36 | c |
37 | c See also line 376 in this file |
38 | call HwU_inithist(nwgt,weights_info) |
39 | - do i=1,2 |
40 | - l=(i-1)*55 |
41 | + do i=1,9 |
42 | + l=(i-1)*59 |
43 | c transverse momenta |
44 | call HwU_book(l+ 1,'total rate '//cc(i),1,0.5d0,1.5d0) |
45 | call HwU_book(l+ 2,'1st charged lepton log pT '//cc(i),25,-0.2d0,3.8d0) |
46 | @@ -76,6 +80,14 @@ |
47 | c HT |
48 | call HwU_book(l+54,'log HT (partons) '//cc(i),25,-0.2d0,3.8d0) |
49 | call HwU_book(l+55,'log HT (reconstructed) '//cc(i),25,-0.2d0,3.8d0) |
50 | + |
51 | + |
52 | + call HwU_book(l+56,'1st isolated ph log pT '//cc(i),25,-0.2d0,3.8d0) |
53 | + call HwU_book(l+57,'2nd isolated ph log pT '//cc(i),25,-0.2d0,3.8d0) |
54 | + call HwU_book(l+58,'3rd isolated ph log pT '//cc(i),25,-0.2d0,3.8d0) |
55 | + |
56 | + call HwU_book(l+59,'3 isolated ph log M '//cc(i),25,-0.2d0,3.8d0) |
57 | + |
58 | enddo |
59 | return |
60 | end |
61 | @@ -110,7 +122,7 @@ |
62 | $ ,ptmpmm,ptepve,ptmmvm,ptepemmpmm,ptepvemmvm,ptt,ptat,pttt |
63 | $ ,etmiss,pth(2),pthh,ptv(3),ptjet(3),Mepem,Mmpmm,Mepve,Mmmvm |
64 | $ ,Mepemmpmm,Mepvemmvm,Mtt,Mhh,Mj1j2,Mj1j3,Mj2j3,Mj1j2j3,Mvvv |
65 | - $ ,HTparton,HTreco,p_reco(0:4,nexternal) |
66 | + $ ,HTparton,HTreco,p_reco(0:4,nexternal),ptphiso(nexternal),Mphphph |
67 | integer nQCD,jet(nexternal),njet,itop,iatop,iem,iep,imp,imm,ive |
68 | $ ,ivm,iv1,iv2,iv3,ih1,ih2,il,ipdg_reco(nexternal) |
69 | double precision getptv4,getptv4_2,getptv4_4,getinvm4_2 |
70 | @@ -119,25 +131,47 @@ |
71 | $ ,getinvm4_4,l10 |
72 | integer orders_tag_plot |
73 | common /corderstagplot/ orders_tag_plot |
74 | + |
75 | + |
76 | +c Photon isolation |
77 | + integer nph,nem,nin,nphiso |
78 | + double precision ptg,chi_gamma_iso,iso_getdrv40 |
79 | + double precision Etsum(0:nexternal) |
80 | + real drlist(nexternal) |
81 | + double precision pgamma(0:3,nexternal),pgammaiso(0:3,nexternal),pem(0:3,nexternal) |
82 | + logical alliso,isolated |
83 | +c Sort array of results: ismode>0 for real, isway=0 for ascending order |
84 | + integer ismode,isway,izero,isorted(nexternal) |
85 | + parameter (ismode=1) |
86 | + parameter (isway=0) |
87 | + parameter (izero=0) |
88 | + integer get_n_tagged_photons |
89 | + |
90 | + |
91 | + logical is_a_lp(nexternal),is_a_lm(nexternal),is_a_j(nexternal) |
92 | + $ ,is_a_ph(nexternal) |
93 | + REAL*8 pt,eta |
94 | + external pt,eta,chi_gamma_iso,sortzv |
95 | +c integer iph1,iph2,iph3 |
96 | c First, try to recombine photons with leptons |
97 | - if (.not.quarkphreco) then |
98 | - write (*,*) 'quark-photon recombination is turned off. '/ |
99 | - $ /'Do need it' |
100 | - stop |
101 | - endif |
102 | - if (.not. lepphreco) then |
103 | - write (*,*) 'lepton-photon recombination is turned off. ' |
104 | - $ //'Do need it.' |
105 | - stop |
106 | - endif |
107 | +c if (.not.quarkphreco) then |
108 | +c write (*,*) 'quark-photon recombination is turned off. '/ |
109 | +c $ /'Do need it' |
110 | +c stop |
111 | +c endif |
112 | +c if (.not. lepphreco) then |
113 | +c write (*,*) 'lepton-photon recombination is turned off. ' |
114 | +c $ //'Do need it.' |
115 | +c stop |
116 | +c endif |
117 | call recombine_momenta(rphreco, etaphreco, lepphreco, quarkphreco, |
118 | $ p, iPDG, p_reco, iPDG_reco) |
119 | |
120 | c Put all (light) QCD partons(+photon) in momentum array for jet clustering. |
121 | nQCD=0 |
122 | do j=nincoming+1,nexternal |
123 | - if (abs(ipdg_reco(j)).le.5 .or. ipdg_reco(j).eq.21 .or. |
124 | - $ ipdg_reco(j).eq.22)then |
125 | + if (abs(ipdg_reco(j)).le.5 .or. ipdg_reco(j).eq.21 |
126 | + $ .or. (ipdg_reco(j).eq.22.and.gamma_is_j)) then |
127 | nQCD=nQCD+1 |
128 | do i=0,3 |
129 | pQCD(i,nQCD)=p_reco(i,j) |
130 | @@ -173,6 +207,161 @@ |
131 | c |
132 | c****************************************************************************** |
133 | |
134 | + |
135 | +c PHOTON (ISOLATION) CUTS |
136 | +c |
137 | +c find the photons |
138 | + do i=1,nexternal |
139 | + if (istatus(i).eq.1 .and. ipdg(i).eq.22 .and. .not.gamma_is_j) then |
140 | + is_a_ph(i)=.true. |
141 | + else |
142 | + is_a_ph(i)=.false. |
143 | + endif |
144 | + enddo |
145 | + if (ptgmin.ne.0d0) then |
146 | + nph=0 |
147 | + do j=nincoming+1,nexternal |
148 | + if (is_a_ph(j)) then |
149 | + nph=nph+1 |
150 | + do i=0,3 |
151 | + pgamma(i,nph)=p(i,j) |
152 | + enddo |
153 | + endif |
154 | + enddo |
155 | + if(nph.eq.0)goto 444 |
156 | +c write(*,*) 'ERROR in cuts.f: photon isolation is not working' |
157 | +c $ // ' for mixed QED-QCD corrections' |
158 | +c stop 1 |
159 | + |
160 | + if(isoEM)then |
161 | + nem=nph |
162 | + do k=1,nem |
163 | + do i=0,3 |
164 | + pem(i,k)=pgamma(i,k) |
165 | + enddo |
166 | + enddo |
167 | + do j=nincoming+1,nexternal |
168 | + if (is_a_lp(j).or.is_a_lm(j)) then |
169 | + nem=nem+1 |
170 | + do i=0,3 |
171 | + pem(i,nem)=p(i,j) |
172 | + enddo |
173 | + endif |
174 | + enddo |
175 | + endif |
176 | + |
177 | +c alliso=.true. |
178 | + nphiso=0 |
179 | + |
180 | + j=0 |
181 | +c do while(j.lt.nph.and.alliso) |
182 | + do while(j.lt.nph) |
183 | + |
184 | +c Loop over all photons |
185 | + j=j+1 |
186 | + |
187 | + ptg=pt(pgamma(0,j)) |
188 | + if(ptg.lt.ptgmin)then |
189 | + cycle |
190 | +c return |
191 | + endif |
192 | + if (etagamma.gt.0d0) then |
193 | + if (abs(eta(pgamma(0,j))).gt.etagamma) then |
194 | + cycle |
195 | +c return |
196 | + endif |
197 | + endif |
198 | + |
199 | +c Isolate from hadronic energy |
200 | + do i=1,nQCD |
201 | + drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pQCD(0,i))) |
202 | + enddo |
203 | + call sortzv(drlist,isorted,nQCD,ismode,isway,izero) |
204 | + Etsum(0)=0.d0 |
205 | + nin=0 |
206 | + do i=1,nQCD |
207 | + if(dble(drlist(isorted(i))).le.R0gamma)then |
208 | + nin=nin+1 |
209 | + Etsum(nin)=Etsum(nin-1)+pt(pQCD(0,isorted(i))) |
210 | + endif |
211 | + enddo |
212 | + isolated=.True. |
213 | + do i=1,nin |
214 | +c alliso=alliso .and. |
215 | +c $ Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))), |
216 | +c $ R0gamma,xn,epsgamma,ptg) |
217 | + if(Etsum(i).gt.chi_gamma_iso(dble(drlist(isorted(i))), |
218 | + $ R0gamma,xn,epsgamma,ptg)) then |
219 | + isolated=isolated.and..False. |
220 | + endif |
221 | + enddo |
222 | + if(.not.isolated)cycle |
223 | + |
224 | +c Isolate from EM energy |
225 | + if(isoEM.and.nem.gt.1)then |
226 | + do i=1,nem |
227 | + drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pem(0,i))) |
228 | + enddo |
229 | + call sortzv(drlist,isorted,nem,ismode,isway,izero) |
230 | +c First of list must be the photon: check this, and drop it |
231 | + if(isorted(1).ne.j.or.drlist(isorted(1)).gt.1.e-4)then |
232 | + write(*,*)'Error #1 in photon isolation' |
233 | + write(*,*)j,isorted(1),drlist(isorted(1)) |
234 | + stop |
235 | + endif |
236 | + Etsum(0)=0.d0 |
237 | + nin=0 |
238 | + do i=2,nem |
239 | + if(dble(drlist(isorted(i))).le.R0gamma)then |
240 | + nin=nin+1 |
241 | + Etsum(nin)=Etsum(nin-1)+pt(pem(0,isorted(i))) |
242 | + endif |
243 | + enddo |
244 | + isolated=.True. |
245 | + do i=1,nin |
246 | +c alliso=alliso .and. |
247 | +c $ Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))), |
248 | +c $ R0gamma,xn,epsgamma,ptg) |
249 | + if(Etsum(i).gt.chi_gamma_iso(dble(drlist(isorted(i))), |
250 | + $ R0gamma,xn,epsgamma,ptg)) then |
251 | + isolated=isolated.and..False. |
252 | + endif |
253 | + enddo |
254 | + if(.not.isolated)cycle |
255 | + endif |
256 | +c End of loop over photons |
257 | + |
258 | + nphiso=nphiso+1 |
259 | + |
260 | + do i=0,3 |
261 | + pgammaiso(i,nphiso)=pgamma(i,j) |
262 | + enddo |
263 | + |
264 | + |
265 | + enddo |
266 | + if(nphiso.lt.get_n_tagged_photons())then |
267 | + print*,"mismatch with cuts.f" |
268 | + stop |
269 | + endif |
270 | + |
271 | + |
272 | +444 continue |
273 | +c End photon isolation |
274 | + endif |
275 | + |
276 | + |
277 | + |
278 | + |
279 | + |
280 | + |
281 | + |
282 | + |
283 | + |
284 | + |
285 | + |
286 | + |
287 | + |
288 | + |
289 | c Look for the other physics objects |
290 | itop=0 |
291 | iatop=0 |
292 | @@ -189,7 +378,14 @@ |
293 | iv3=0 |
294 | ih1=0 |
295 | ih2=0 |
296 | +c iph1=0 |
297 | +c iph2=0 |
298 | +c iph3=0 |
299 | +c print*,"nell'analisi" |
300 | do i=1,nexternal |
301 | +c print*,"idpg di ",i,"=",ipdg(i) |
302 | +c print*,"idpg_reco di ",i,"=",ipdg_reco(i) |
303 | + |
304 | if (ipdg_reco(i).eq.6) then |
305 | itop=i |
306 | elseif(ipdg_reco(i).eq.-6) then |
307 | @@ -232,8 +428,27 @@ |
308 | stop |
309 | endif |
310 | endif |
311 | +c elseif(abs(ipdg_reco(i)).eq.22.and.i.gt.nincoming) then |
312 | +c if (iph1.eq.0) then |
313 | +c iph1=i |
314 | +c else |
315 | +c if (iph2.eq.0) then |
316 | +c iph2=i |
317 | +c else |
318 | +c if (iph3.eq.0) then |
319 | +c iph3=i |
320 | +c else |
321 | +c write (*,*) 'too many photonss' |
322 | +c stop |
323 | +c endif |
324 | +c endif |
325 | +c endif |
326 | endif |
327 | + |
328 | enddo |
329 | +c print*,itop,iatop |
330 | +c stop |
331 | + |
332 | if (itop.ne.0) ptt=getptv4(p_reco(0,itop)) |
333 | if (iatop.ne.0) ptat=getptv4(p_reco(0,iatop)) |
334 | if (itop.ne.0 .and. iatop.ne.0) then |
335 | @@ -318,11 +533,21 @@ |
336 | pth(2)=tmp |
337 | endif |
338 | endif |
339 | + |
340 | + |
341 | if (iv1.ne.0) ptv(1)=getptv4(p_reco(0,iv1)) |
342 | if (iv2.ne.0) ptv(2)=getptv4(p_reco(0,iv2)) |
343 | if (iv3.ne.0) ptv(3)=getptv4(p_reco(0,iv3)) |
344 | if (iv1.ne.0 .and. iv2.ne.0 .and. iv3.ne.0) then |
345 | Mvvv=getinvm4_3(p_reco(0,iv1),p_reco(0,iv2),p_reco(0,iv3)) |
346 | + endif |
347 | + |
348 | + do i=1,nphiso |
349 | + ptphiso(i)=getptv4(pgammaiso(0,i)) |
350 | + enddo |
351 | + |
352 | + if (iv1.ne.0 .and. iv2.ne.0 .and. iv3.ne.0) then |
353 | + Mvvv=getinvm4_3(p_reco(0,iv1),p_reco(0,iv2),p_reco(0,iv3)) |
354 | c order the vector bosons (if there are 3) |
355 | do i=1,2 |
356 | do j=1,3-i |
357 | @@ -341,6 +566,31 @@ |
358 | ptv(2)=tmp |
359 | endif |
360 | endif |
361 | + |
362 | + |
363 | + |
364 | + |
365 | + if (nphiso.eq.3) then |
366 | + Mphphph=getinvm4_3(pgammaiso(0,1),pgammaiso(0,2),pgammaiso(0,3)) |
367 | +c order the isolated photons (if there are 3) |
368 | + do i=1,2 |
369 | + do j=1,3-i |
370 | + if (ptphiso(j).lt.ptphiso(j+1)) then |
371 | + tmp=ptphiso(j) |
372 | + ptphiso(j)=ptphiso(j+1) |
373 | + ptphiso(j+1)=tmp |
374 | + endif |
375 | + enddo |
376 | + enddo |
377 | + elseif (nphiso.eq.2) then |
378 | +c order the isolated photons (if there are 2) |
379 | + if (ptphiso(1).lt.ptphiso(2)) then |
380 | + tmp=ptphiso(1) |
381 | + ptphiso(1)=ptphiso(2) |
382 | + ptphiso(2)=tmp |
383 | + endif |
384 | + endif |
385 | + |
386 | do i=1,njet |
387 | ptjet(i)=getptv4(pjet(0,i)) |
388 | enddo |
389 | @@ -369,9 +619,18 @@ |
390 | enddo |
391 | if (ive.ne.0 .or. ivm.ne.0) HTreco=HTreco+etmiss |
392 | |
393 | - do i=1,2 |
394 | - l=(i-1)*55 |
395 | + do i=1,9 |
396 | + l=(i-1)*59 |
397 | if (ibody.ne.3 .and.i.eq.2) cycle |
398 | + if (i.eq. 3.and.orders_tag_plot.ne.204) cycle |
399 | + if (i.eq. 4.and.orders_tag_plot.ne.402) cycle |
400 | + if (i.eq. 5.and.orders_tag_plot.ne.600) cycle |
401 | + if (i.eq. 6.and.orders_tag_plot.ne.206) cycle |
402 | + if (i.eq. 7.and.orders_tag_plot.ne.404) cycle |
403 | + if (i.eq. 8.and.orders_tag_plot.ne.602) cycle |
404 | + if (i.eq. 9.and.orders_tag_plot.ne.800) cycle |
405 | + |
406 | + |
407 | |
408 | c How to tag orders (QCD+QED*100) |
409 | c |
410 | @@ -493,6 +752,18 @@ |
411 | c HT |
412 | call HwU_fill(l+54,l10(HTparton),wgts) |
413 | call HwU_fill(l+55,l10(HTreco),wgts) |
414 | + |
415 | + if (nphiso.ge.1) call HwU_fill(l+56,l10(ptphiso(1)),wgts) |
416 | + if (nphiso.ge.2) call HwU_fill(l+57,l10(ptphiso(2)),wgts) |
417 | + if (nphiso.ge.3) call HwU_fill(l+58,l10(ptphiso(3)),wgts) |
418 | + |
419 | + if (nphiso.ge.3) call HwU_fill(l+59,l10(Mphphph),wgts) |
420 | + |
421 | + |
422 | + |
423 | + |
424 | + |
425 | + |
426 | enddo |
427 | |
428 | 999 return |
429 | |
430 | === modified file 'Template/NLO/SubProcesses/BinothLHA.f' |
431 | --- Template/NLO/SubProcesses/BinothLHA.f 2020-06-26 09:54:58 +0000 |
432 | +++ Template/NLO/SubProcesses/BinothLHA.f 2021-10-15 11:42:16 +0000 |
433 | @@ -8,6 +8,7 @@ |
434 | implicit none |
435 | include "nexternal.inc" |
436 | include "coupl.inc" |
437 | + include "../../Source/MODEL/input.inc" |
438 | include 'born_nhel.inc' |
439 | double precision pi, zero,mone |
440 | parameter (pi=3.1415926535897932385d0) |
441 | @@ -57,6 +58,7 @@ |
442 | integer getordpowfromindex_ml5 |
443 | integer orders_to_amp_split_pos |
444 | logical, allocatable, save :: keep_order(:) |
445 | + include 'a0Gmuconv.inc' |
446 | include 'orders.inc' |
447 | integer amp_orders(nsplitorders) |
448 | integer split_amp_orders(nsplitorders), iamp |
449 | @@ -74,10 +76,17 @@ |
450 | |
451 | logical updateloop |
452 | common /to_updateloop/updateloop |
453 | + |
454 | + integer get_n_tagged_photons |
455 | + external get_n_tagged_photons |
456 | + double precision get_virtual_a0Gmu_conv |
457 | + external get_virtual_a0Gmu_conv |
458 | + integer ntagph, qed_pow_b |
459 | c masses |
460 | include 'pmass.inc' |
461 | data nbad / 0 / |
462 | |
463 | + |
464 | IOErrCounter = 0 |
465 | c update the ren_scale for MadLoop and the couplings (should be the |
466 | c Ellis-Sexton scale) |
467 | @@ -93,8 +102,7 @@ |
468 | single = 0d0 |
469 | double = 0d0 |
470 | born_hel_from_virt = 0d0 |
471 | -C reset the amp_split array |
472 | - amp_split(1:amp_split_size) = 0d0 |
473 | +C reset the various arrays |
474 | amp_split_finite_ML(1:amp_split_size) = 0d0 |
475 | amp_split_poles_ML(1:amp_split_size,1) = 0d0 |
476 | amp_split_poles_ML(1:amp_split_size,2) = 0d0 |
477 | @@ -256,6 +264,32 @@ |
478 | c endif |
479 | c virt_wgt=virt_wgt+conversion*born_wgt*ao2pi |
480 | c====================================================================== |
481 | + |
482 | +c====================================================================== |
483 | +c If there are tagged photon and other particles in the process, |
484 | +C one must use a mixed Gmu-alpha0 renormalisation. |
485 | + ntagph = get_n_tagged_photons() |
486 | + if (ntagph.ne.0) then |
487 | + do i = 1, AMP_SPLIT_SIZE_BORN |
488 | + call amp_split_pos_to_orders(i, amp_orders) |
489 | + born_wgt = amp_split(i) |
490 | + ! this is the number of powers of 'e' in the born |
491 | + qed_pow_b = amp_orders(qed_pos) |
492 | + amp_orders(qed_pos) = amp_orders(qed_pos) + 2 |
493 | + if (amp_orders(qed_pos).gt.NLO_ORDERS(qed_pos)) cycle |
494 | + |
495 | + amp_split_poles_ML(orders_to_amp_split_pos(amp_orders),1) = |
496 | + $ amp_split_poles_ML(orders_to_amp_split_pos(amp_orders),1) + |
497 | + $ get_virtual_a0Gmu_conv(qed_pow_b,ntagph,1,born_wgt) |
498 | + |
499 | + amp_split_finite_ML(orders_to_amp_split_pos(amp_orders)) = |
500 | + $ amp_split_finite_ML(orders_to_amp_split_pos(amp_orders)) + |
501 | + $ get_virtual_a0Gmu_conv(qed_pow_b,ntagph,0,born_wgt) |
502 | + |
503 | + enddo |
504 | + endif |
505 | +c====================================================================== |
506 | + |
507 | c |
508 | c Check poles for the first PS points when doing MC over helicities, and |
509 | c for all phase-space points when not doing MC over helicities. Skip |
510 | @@ -268,7 +302,7 @@ |
511 | polecheck_passed = .true. |
512 | ! loop over the full result and each of the amp_split |
513 | ! contribution |
514 | - do iamp=0,amp_split_size |
515 | + do iamp=1,amp_split_size |
516 | ! skip 0 contributions in the amp_split array |
517 | if (iamp.ne.0) then |
518 | if (amp_split_poles_FKS(iamp,1).eq.0d0.and. |
519 | |
520 | === modified file 'Template/NLO/SubProcesses/chooser_functions.f' |
521 | --- Template/NLO/SubProcesses/chooser_functions.f 2018-04-05 17:59:21 +0000 |
522 | +++ Template/NLO/SubProcesses/chooser_functions.f 2021-10-15 11:42:16 +0000 |
523 | @@ -100,6 +100,10 @@ |
524 | double precision ch_i,ch_j,ch_m |
525 | integer particle_type_born(nexternal-1) |
526 | common /c_particle_type_born/particle_type_born |
527 | + logical particle_tag(nexternal) |
528 | + common /c_particle_tag/particle_tag |
529 | + logical particle_tag_born(nexternal-1) |
530 | + common /c_particle_tag/particle_tag_born |
531 | logical need_color_links, need_charge_links |
532 | common /c_need_links/need_color_links, need_charge_links |
533 | integer extra_cnt, isplitorder_born, isplitorder_cnt |
534 | @@ -130,6 +134,7 @@ |
535 | stop |
536 | endif |
537 | particle_type(i)=particle_type_D(nFKSprocess,i) |
538 | + particle_tag(i)=particle_tag_D(nFKSprocess,i) |
539 | particle_charge(i)=particle_charge_D(nFKSprocess,i) |
540 | pdg_type(i)=pdg_type_D(nFKSprocess,i) |
541 | ! is_aorg is true if the particle can induce soft singularities |
542 | @@ -144,18 +149,22 @@ |
543 | if (i.lt.min(i_fks,j_fks)) then |
544 | particle_type_born(i)=particle_type(i) |
545 | particle_charge_born(i)=particle_charge(i) |
546 | + particle_tag_born(i)=particle_tag(i) |
547 | elseif (i.gt.max(i_fks,j_fks)) then |
548 | particle_type_born(i-1)=particle_type(i) |
549 | particle_charge_born(i-1)=particle_charge(i) |
550 | + particle_tag_born(i-1)=particle_tag(i) |
551 | elseif (i.eq.min(i_fks,j_fks)) then |
552 | i_type=particle_type(i_fks) |
553 | j_type=particle_type(j_fks) |
554 | ch_i=particle_charge(i_fks) |
555 | ch_j=particle_charge(j_fks) |
556 | + particle_tag_born(i) = particle_tag(j_fks) |
557 | call get_mother_col_charge(i_type,ch_i,j_type,ch_j,m_type,ch_m) |
558 | particle_type_born(i)=m_type |
559 | particle_charge_born(i)=ch_m |
560 | elseif (i.ne.max(i_fks,j_fks)) then |
561 | + particle_tag_born(i) = particle_tag(i) |
562 | particle_type_born(i)=particle_type(i) |
563 | particle_charge_born(i)=particle_charge(i) |
564 | endif |
565 | |
566 | === modified file 'Template/NLO/SubProcesses/cuts.f' |
567 | --- Template/NLO/SubProcesses/cuts.f 2021-01-08 10:23:30 +0000 |
568 | +++ Template/NLO/SubProcesses/cuts.f 2021-10-15 11:42:16 +0000 |
569 | @@ -40,42 +40,28 @@ |
570 | external R2_04,invm2_04,pt_04,eta_04,pt,eta |
571 | C recombination of photons |
572 | double precision p_reco(0:4,nexternal), R_reco |
573 | - integer iPDG_reco(nexternal) |
574 | + integer iPDG_reco(nexternal),nphiso |
575 | c local integers |
576 | integer i,j |
577 | -c temporary variable for caching locally computation |
578 | - double precision tmpvar |
579 | +c bare parton algorithm |
580 | + integer nPART |
581 | + double precision pPART(0:3,nexternal) |
582 | c jet cluster algorithm |
583 | integer nQCD,NJET,JET(nexternal) |
584 | double precision pQCD(0:3,nexternal),PJET(0:3,nexternal) |
585 | - double precision rfj,sycut,palg,amcatnlo_fastjetdmerge |
586 | integer njet_eta |
587 | - integer mm |
588 | -c Photon isolation |
589 | - integer nph,nem,k,nin |
590 | - double precision ptg,chi_gamma_iso,iso_getdrv40 |
591 | - double precision Etsum(0:nexternal) |
592 | - real drlist(nexternal) |
593 | - double precision pgamma(0:3,nexternal),pem(0:3,nexternal) |
594 | - logical alliso |
595 | -c Sort array of results: ismode>0 for real, isway=0 for ascending order |
596 | - integer ismode,isway,izero,isorted(nexternal) |
597 | - parameter (ismode=1) |
598 | - parameter (isway=0) |
599 | - parameter (izero=0) |
600 | -c The UNLOPS cut |
601 | - double precision p_unlops(0:3,nexternal) |
602 | + double precision pgamma(0:3,nexternal),pgamma_iso(0:3,nexternal) |
603 | + integer nph |
604 | include "run.inc" ! includes the ickkw parameter |
605 | - logical passUNLOPScuts |
606 | -c PDG specific cut |
607 | - double precision etmin(nincoming+1:nexternal-1) |
608 | - double precision etmax(nincoming+1:nexternal-1) |
609 | - double precision mxxmin(nincoming+1:nexternal-1,nincoming+1:nexternal-1) |
610 | - common /to_cuts/etmin,etmax,mxxmin |
611 | c logicals that define if particles are leptons, jets or photons. These |
612 | c are filled from the PDG codes (iPDG array) in this function. |
613 | logical is_a_lp(nexternal),is_a_lm(nexternal),is_a_j(nexternal) |
614 | - $ ,is_a_ph(nexternal) |
615 | + $,is_a_ph(nexternal),is_nph_iso(nexternal),is_nextph_iso(nexternal) |
616 | + $,is_nextph_iso_reco(nexternal) |
617 | + logical is_a_lp_reco(nexternal),is_a_lm_reco(nexternal) |
618 | + logical passcuts_leptons, passcuts_unlops_jv, passcuts_photons, |
619 | + $ passcuts_jets, passcuts_pdgs |
620 | + |
621 | |
622 | passcuts_user=.true. ! event is okay; otherwise it is changed |
623 | |
624 | @@ -84,98 +70,329 @@ |
625 | C Cuts from the run_card.dat |
626 | C*************************************************************** |
627 | C*************************************************************** |
628 | - !first recombine the photons and fermions |
629 | + |
630 | + ! Find the bare QCD partons |
631 | + ! This is used only as input to the photon isolation |
632 | + call identify_PART_partons(p,istatus,ipdg,pPART,nPART,is_a_lp,is_a_lm) |
633 | + |
634 | + ! Apply the Photon cuts on isolated photons based on the bare particles |
635 | + passcuts_user = passcuts_user .and. |
636 | + $ passcuts_photons(p,istatus,ipdg,is_a_lp,is_a_lm,pPART,nPART, |
637 | + $ pgamma,nph,is_nph_iso,is_nextph_iso) |
638 | + if (.not.passcuts_user) return |
639 | + |
640 | + ! Recombine the photons and fermions |
641 | call recombine_momenta(rphreco, etaphreco, lepphreco, quarkphreco, |
642 | - $ p, iPDG, p_reco, iPDG_reco) |
643 | - |
644 | -c |
645 | -c CHARGED LEPTON CUTS |
646 | -c |
647 | -c find the charged leptons (also used in the photon isolation cuts below) |
648 | - do i=1,nexternal |
649 | - if(istatus(i).eq.1 .and. |
650 | - & (ipdg_reco(i).eq.11 .or. ipdg_reco(i).eq.13 .or. ipdg_reco(i).eq.15)) then |
651 | - is_a_lm(i)=.true. |
652 | - else |
653 | - is_a_lm(i)=.false. |
654 | - endif |
655 | - if(istatus(i).eq.1 .and. |
656 | - & (ipdg_reco(i).eq.-11 .or. ipdg_reco(i).eq.-13 .or. ipdg_reco(i).eq.-15)) then |
657 | - is_a_lp(i)=.true. |
658 | - else |
659 | - is_a_lp(i)=.false. |
660 | - endif |
661 | - enddo |
662 | -c apply the charged lepton cuts |
663 | - do i=nincoming+1,nexternal |
664 | - if (is_a_lp(i).or.is_a_lm(i)) then |
665 | -c transverse momentum |
666 | - if (ptl.gt.0d0) then |
667 | - if (pt_04(p_reco(0,i)).lt.ptl) then |
668 | - passcuts_user=.false. |
669 | - return |
670 | - endif |
671 | - endif |
672 | -c pseudo-rapidity |
673 | - if (etal.gt.0d0) then |
674 | - if (abs(eta_04(p_reco(0,i))).gt.etal) then |
675 | - passcuts_user=.false. |
676 | - return |
677 | - endif |
678 | - endif |
679 | -c DeltaR and invariant mass cuts |
680 | - if (is_a_lp(i)) then |
681 | - do j=nincoming+1,nexternal |
682 | - if (is_a_lm(j)) then |
683 | - if (drll.gt.0d0) then |
684 | - if (R2_04(p_reco(0,i),p_reco(0,j)).lt.drll**2) then |
685 | - passcuts_user=.false. |
686 | - return |
687 | - endif |
688 | - endif |
689 | - if (mll.gt.0d0) then |
690 | - if (invm2_04(p_reco(0,i),p_reco(0,j),1d0).lt.mll**2) then |
691 | - passcuts_user=.false. |
692 | - return |
693 | - endif |
694 | - endif |
695 | - if (ipdg_reco(i).eq.-ipdg_reco(j)) then |
696 | - if (drll_sf.gt.0d0) then |
697 | - if (R2_04(p_reco(0,i),p_reco(0,j)).lt.drll_sf**2) then |
698 | - passcuts_user=.false. |
699 | - return |
700 | - endif |
701 | - endif |
702 | - if (mll_sf.gt.0d0) then |
703 | - if (invm2_04(p_reco(0,i),p_reco(0,j),1d0).lt.mll_sf**2) |
704 | - $ then |
705 | - passcuts_user=.false. |
706 | - return |
707 | - endif |
708 | - endif |
709 | - endif |
710 | + $ p, iPDG, is_nextph_iso, p_reco, iPDG_reco, is_nextph_iso_reco) |
711 | + |
712 | + ! Apply the reco lepton cuts |
713 | + passcuts_user = passcuts_user .and. |
714 | + $ passcuts_leptons(p_reco,istatus,ipdg_reco,is_a_lp_reco,is_a_lm_reco) |
715 | + if (.not.passcuts_user) return |
716 | + |
717 | + ! Find the reco QCD partons including |
718 | + ! A. All photons if gamma_is_j is on |
719 | + ! B. Non-iso, non-reco photons if reco is on |
720 | + call identify_QCD_partons(is_nextph_iso_reco,p_reco,istatus,ipdg_reco,is_a_j,pQCD,nQCD) |
721 | + |
722 | + ! Apply the UNLOPS/JetVeto cuts |
723 | + passcuts_user = passcuts_user .and. |
724 | + $ passcuts_unlops_jv(p_reco,istatus,ipdg_reco,pQCD,nQCD,ickkw) |
725 | + if (.not.passcuts_user) return |
726 | + |
727 | + ! Apply the Jet cuts |
728 | + passcuts_user = passcuts_user .and. |
729 | + $ passcuts_jets(p_reco,pQCD,nQCD,pgamma,nph,is_nph_iso,ickkw) |
730 | + if (.not.passcuts_user) return |
731 | + |
732 | + ! Apply PDG specific cuts |
733 | + passcuts_user = passcuts_user .and. |
734 | + $ passcuts_pdgs(p_reco,istatus,ipdg_reco) |
735 | + if (.not.passcuts_user) return |
736 | + |
737 | +C*************************************************************** |
738 | +C*************************************************************** |
739 | +C PUT HERE YOUR USER-DEFINED CUTS |
740 | +C*************************************************************** |
741 | +C*************************************************************** |
742 | +C |
743 | +c$$$C EXAMPLE: cut on top quark pT |
744 | +c$$$C Note that PDG specific cut are more optimised than simple user cut |
745 | +c$$$ do i=1,nexternal ! loop over all external particles |
746 | +c$$$ if (istatus(i).eq.1 ! final state particle |
747 | +c$$$ & .and. abs(ipdg(i)).eq.6) then ! top quark |
748 | +c$$$C apply the pT cut (pT should be large than 200 GeV for the event to |
749 | +c$$$C pass cuts) |
750 | +c$$$ if ( p(1,i)**2+p(2,i)**2 .lt. 200d0**2 ) then |
751 | +c$$$C momenta do not pass cuts. Set passcuts_user to false and return |
752 | +c$$$ passcuts_user=.false. |
753 | +c$$$ return |
754 | +c$$$ endif |
755 | +c$$$ endif |
756 | +c$$$ enddo |
757 | +c |
758 | + return |
759 | + end |
760 | + |
761 | + subroutine identify_PART_partons(p,istatus,ipdg,pPART,nPART,is_a_lp,is_a_lm) |
762 | + implicit none |
763 | + include 'nexternal.inc' |
764 | + integer istatus(nexternal) |
765 | + integer iPDG(nexternal) |
766 | + double precision p(0:4,nexternal) |
767 | + integer nPART |
768 | + double precision pPART(0:3,nexternal) |
769 | + logical is_a_lp(nexternal),is_a_lm(nexternal) |
770 | + include "run.inc" |
771 | + include "cuts.inc" |
772 | + |
773 | + integer i, j |
774 | +c |
775 | +c Bare partons and leptons |
776 | +c |
777 | + nPART=0 |
778 | + do j=1,nexternal |
779 | + is_a_lp(j)=.false. |
780 | + is_a_lm(j)=.false. |
781 | +c Partons |
782 | + if (istatus(j).eq.1 .and. |
783 | + & (abs(ipdg(j)).le.maxjetflavor .or. ipdg(j).eq.21) |
784 | + &) then |
785 | + nPART=nPART+1 |
786 | + do i=0,3 |
787 | + pPART(i,nPART)=p(i,j) |
788 | + enddo |
789 | + endif |
790 | +c Leptons |
791 | + if (ipdg(j).eq.11 .or. ipdg(j).eq.13 |
792 | + $ .or. ipdg(j).eq.15) then |
793 | + is_a_lm(j)=.true. |
794 | + endif |
795 | + if (ipdg(j).eq.-11 .or. ipdg(j).eq.-13 |
796 | + $ .or. ipdg(j).eq.-15) then |
797 | + is_a_lp(j)=.true. |
798 | + endif |
799 | + enddo |
800 | + |
801 | + return |
802 | + end |
803 | + |
804 | + logical function passcuts_photons(p,istatus,ipdg,is_a_lp,is_a_lm, |
805 | + $pPART,nPART,pgamma,nph,is_nph_iso,is_nextph_iso) |
806 | + implicit none |
807 | + include 'nexternal.inc' |
808 | + integer istatus(nexternal) |
809 | + integer iPDG(nexternal) |
810 | + double precision p(0:4,nexternal) |
811 | + logical is_a_lp(nexternal),is_a_lm(nexternal) |
812 | + integer nPART, nph |
813 | + double precision pPART(0:3,nexternal), pgamma(0:3,nexternal) |
814 | + double precision pgamma_iso(0:3,nexternal) |
815 | + logical is_nph_iso(nexternal),is_nextph_iso(nexternal) |
816 | + include "cuts.inc" |
817 | + include "run.inc" |
818 | + integer i,j,k,mu |
819 | +c Sort array of results: ismode>0 for real, isway=0 for ascending order |
820 | + integer ismode,isway,izero,isorted(nexternal) |
821 | + parameter (ismode=1) |
822 | + parameter (isway=0) |
823 | + parameter (izero=0) |
824 | + |
825 | +c Photon isolation |
826 | + integer nem,nin,nphiso |
827 | + double precision ptg,chi_gamma_iso,iso_getdrv40 |
828 | + double precision Etsum(0:nexternal) |
829 | + real drlist(nexternal) |
830 | + double precision pem(0:3,nexternal) |
831 | + |
832 | + logical alliso,isolated |
833 | + integer get_n_tagged_photons |
834 | + logical is_a_ph(nexternal) |
835 | + |
836 | + REAL*8 pt,eta |
837 | + external pt,eta |
838 | + |
839 | + passcuts_photons = .true. |
840 | + |
841 | +c |
842 | +c PHOTON (ISOLATION) CUTS |
843 | +c |
844 | +c Initialise common logical iso |
845 | + do i=nincoming+1,nexternal |
846 | + is_nextph_iso(i)=.False. |
847 | + enddo |
848 | +c find the photons |
849 | + do i=nincoming+1,nexternal |
850 | + if (ipdg(i).eq.22 .and. .not.gamma_is_j) then |
851 | + is_a_ph(i)=.true. |
852 | + else |
853 | + is_a_ph(i)=.false. |
854 | + endif |
855 | + enddo |
856 | + |
857 | + if (ptgmin.ne.0d0) then |
858 | + nph=0 |
859 | + do j=nincoming+1,nexternal |
860 | + if (is_a_ph(j)) then |
861 | + nph=nph+1 |
862 | + do i=0,3 |
863 | + pgamma(i,nph)=p(i,j) |
864 | + enddo |
865 | + endif |
866 | + enddo |
867 | + if(nph.eq.0) return |
868 | + |
869 | + if(isoEM)then |
870 | + nem=nph |
871 | + do k=1,nem |
872 | + do i=0,3 |
873 | + pem(i,k)=pgamma(i,k) |
874 | + enddo |
875 | + enddo |
876 | + do j=nincoming+1,nexternal |
877 | + if (is_a_lp(j).or.is_a_lm(j)) then |
878 | + nem=nem+1 |
879 | + do i=0,3 |
880 | + pem(i,nem)=p(i,j) |
881 | + enddo |
882 | + endif |
883 | + enddo |
884 | + endif |
885 | + |
886 | + nphiso=0 |
887 | + |
888 | + j=0 |
889 | +c Loop over all photons |
890 | + do while(j.lt.nph) |
891 | + |
892 | + j=j+1 |
893 | + is_nph_iso(j)=.False. |
894 | + ptg=pt(pgamma(0,j)) |
895 | + if(ptg.lt.ptgmin)then |
896 | + cycle |
897 | + endif |
898 | + if (etagamma.gt.0d0) then |
899 | + if (abs(eta(pgamma(0,j))).gt.etagamma) then |
900 | + cycle |
901 | + endif |
902 | + endif |
903 | + |
904 | +c Isolate from hadronic energy |
905 | + do i=1,nPART |
906 | + drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pPART(0,i))) |
907 | + enddo |
908 | + call sortzv(drlist,isorted,nPART,ismode,isway,izero) |
909 | + Etsum(0)=0.d0 |
910 | + nin=0 |
911 | + do i=1,nPART |
912 | + if(dble(drlist(isorted(i))).le.R0gamma)then |
913 | + nin=nin+1 |
914 | + Etsum(nin)=Etsum(nin-1)+pt(pPART(0,isorted(i))) |
915 | + endif |
916 | + enddo |
917 | + isolated=.True. |
918 | + do i=1,nin |
919 | + if(Etsum(i).gt.chi_gamma_iso(dble(drlist(isorted(i))), |
920 | + $ R0gamma,xn,epsgamma,ptg)) then |
921 | + isolated=.False. |
922 | + exit |
923 | + endif |
924 | + enddo |
925 | + if(.not.isolated)cycle |
926 | + |
927 | +c Isolate from EM energy |
928 | + if(isoEM.and.nem.gt.1)then |
929 | + do i=1,nem |
930 | + drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pem(0,i))) |
931 | + enddo |
932 | + call sortzv(drlist,isorted,nem,ismode,isway,izero) |
933 | +c First of list must be the photon: check this, and drop it |
934 | + if(isorted(1).ne.j.or.drlist(isorted(1)).gt.1.e-4)then |
935 | + write(*,*)'Error #1 in photon isolation' |
936 | + write(*,*)j,isorted(1),drlist(isorted(1)) |
937 | + stop |
938 | + endif |
939 | + Etsum(0)=0.d0 |
940 | + nin=0 |
941 | + do i=2,nem |
942 | + if(dble(drlist(isorted(i))).le.R0gamma)then |
943 | + nin=nin+1 |
944 | + Etsum(nin)=Etsum(nin-1)+pt(pem(0,isorted(i))) |
945 | endif |
946 | enddo |
947 | + isolated=.True. |
948 | + do i=1,nin |
949 | + if(Etsum(i).gt.chi_gamma_iso(dble(drlist(isorted(i))), |
950 | + $ R0gamma,xn,epsgamma,ptg)) then |
951 | + isolated=.False. |
952 | + exit |
953 | + endif |
954 | + enddo |
955 | + if(.not.isolated)cycle |
956 | endif |
957 | + is_nph_iso(j)=.True. |
958 | + nphiso=nphiso+1 |
959 | + |
960 | + if (nphiso.gt.0) then |
961 | + do mu=0,3 |
962 | + pgamma_iso(mu,nphiso)=pgamma(mu,j) |
963 | + enddo |
964 | + |
965 | + do i=nincoming+1,nexternal |
966 | + if ( ipdg(i).eq.22 .and. |
967 | + $ pt(p(0,i)).eq.pt(pgamma_iso(0,nphiso)) ) then |
968 | + is_nextph_iso(i)=.True. |
969 | + endif |
970 | + enddo |
971 | + endif |
972 | + enddo |
973 | +c End of loop over photons |
974 | + |
975 | + if(nphiso.lt.get_n_tagged_photons())then |
976 | + passcuts_photons=.false. |
977 | + return |
978 | endif |
979 | - enddo |
980 | + endif |
981 | + |
982 | + return |
983 | + end |
984 | + |
985 | + |
986 | + subroutine identify_QCD_partons(is_iso,p,istatus,ipdg,is_a_j,pQCD,nQCD) |
987 | + implicit none |
988 | + include 'nexternal.inc' |
989 | + integer istatus(nexternal) |
990 | + integer iPDG(nexternal) |
991 | + double precision p(0:4,nexternal) |
992 | + logical is_a_j(nexternal) |
993 | + integer nQCD |
994 | + double precision pQCD(0:3,nexternal) |
995 | + logical is_iso(nexternal) |
996 | + REAL*8 pt,eta |
997 | + external pt,eta |
998 | + include "run.inc" |
999 | + include "cuts.inc" |
1000 | + |
1001 | + integer i, j |
1002 | + |
1003 | c |
1004 | c JET CUTS |
1005 | c |
1006 | c find the jets |
1007 | do i=1,nexternal |
1008 | if (istatus(i).eq.1 .and. |
1009 | - & (abs(ipdg_reco(i)).le.maxjetflavor .or. ipdg_reco(i).eq.21 |
1010 | - & .or.(ipdg_reco(i).eq.22.and.gamma_is_j))) then |
1011 | + & ( abs(ipdg(i)).le.maxjetflavor .or. ipdg(i).eq.21 |
1012 | + & .or. (ipdg(i).eq.22.and.gamma_is_j) .or. |
1013 | + & (ipdg(i).eq.22.and. .not.is_iso(i)) ) |
1014 | + &) then |
1015 | is_a_j(i)=.true. |
1016 | else |
1017 | is_a_j(i)=.false. |
1018 | endif |
1019 | enddo |
1020 | - |
1021 | c If we do not require a mimimum jet energy, there's no need to apply |
1022 | c jet clustering and all that. |
1023 | - if (ptj.gt.0d0.or.ptgmin.gt.0d0) then |
1024 | + if (ptj.ne.0d0.or.ptgmin.ne.0d0) then |
1025 | c Put all (light) QCD partons in momentum array for jet clustering. |
1026 | c From the run_card.dat, maxjetflavor defines if b quark should be |
1027 | c considered here (via the logical variable 'is_a_jet'). nQCD becomes |
1028 | @@ -186,41 +403,43 @@ |
1029 | if (is_a_j(j)) then |
1030 | nQCD=nQCD+1 |
1031 | do i=0,3 |
1032 | - pQCD(i,nQCD)=p_reco(i,j) |
1033 | + pQCD(i,nQCD)=p(i,j) |
1034 | enddo |
1035 | endif |
1036 | enddo |
1037 | endif |
1038 | |
1039 | -c THE UNLOPS CUT: |
1040 | - if (ickkw.eq.4 .and. ptj.gt.0d0) then |
1041 | -c Use special pythia pt cut for minimal pT |
1042 | - do i=1,nexternal |
1043 | - do j=0,3 |
1044 | - p_unlops(j,i)=p_reco(j,i) |
1045 | - enddo |
1046 | - enddo |
1047 | - call pythia_UNLOPS(p_unlops,passUNLOPScuts) |
1048 | - if (.not. passUNLOPScuts) then |
1049 | - passcuts_user=.false. |
1050 | - return |
1051 | - endif |
1052 | -c Bypass normal jet cuts |
1053 | - goto 122 |
1054 | -c THE VETO XSEC CUT: |
1055 | - elseif (ickkw.eq.-1 .and. ptj.gt.0d0) then |
1056 | -c Use veto'ed Xsec for analytic NNLL resummation |
1057 | - if (nQCD.ne.1) then |
1058 | - write (*,*) 'ERROR: more than one QCD parton in '/ |
1059 | - $ /'this event in cuts.f. There should only be one' |
1060 | - stop |
1061 | - endif |
1062 | - if (pt(pQCD(0,1)) .gt. ptj) then |
1063 | - passcuts_user=.false. |
1064 | - return |
1065 | - endif |
1066 | - endif |
1067 | - |
1068 | + return |
1069 | + end |
1070 | + |
1071 | + |
1072 | + logical function passcuts_jets(p,pQCD,nQCD,pgamma,nph,is_nph_iso,ickkw) |
1073 | + implicit none |
1074 | + include 'nexternal.inc' |
1075 | + double precision p(0:4,nexternal) |
1076 | + integer nQCD, nph |
1077 | + double precision pQCD(0:3,nexternal), pgamma(0:3,nexternal) |
1078 | + logical is_nph_iso(nexternal) |
1079 | + integer ickkw |
1080 | + include "cuts.inc" |
1081 | + |
1082 | + integer NJET,JET(nexternal) |
1083 | + double precision rfj,sycut,palg,amcatnlo_fastjetdmerge |
1084 | + double precision PJET(0:3,nexternal) |
1085 | + integer mm |
1086 | + integer i,j |
1087 | + |
1088 | + integer get_n_tagged_photons |
1089 | + |
1090 | + REAL*8 pt,eta |
1091 | + external pt,eta |
1092 | + |
1093 | + passcuts_jets=.true. |
1094 | + |
1095 | +c JET CUTS |
1096 | + |
1097 | +C do nothing if ickkw=4 (UNLOPS) |
1098 | + if (ickkw.eq.4)return |
1099 | |
1100 | if (ptj.gt.0d0.and.nQCD.gt.1) then |
1101 | |
1102 | @@ -233,11 +452,10 @@ |
1103 | if(abs(pQCD(0,j)/p(0,1)).lt.1.d-8) mm=mm+1 |
1104 | enddo |
1105 | if(mm.gt.1)then |
1106 | - passcuts_user=.false. |
1107 | + passcuts_jets=.false. |
1108 | return |
1109 | endif |
1110 | |
1111 | - |
1112 | c Define jet clustering parameters (from cuts.inc via the run_card.dat) |
1113 | palg=JETALGO ! jet algorithm: 1.0=kt, 0.0=C/A, -1.0 = anti-kt |
1114 | rfj=JETRADIUS ! the radius parameter |
1115 | @@ -268,127 +486,178 @@ |
1116 | |
1117 | c Apply the jet cuts |
1118 | if (njet .ne. nQCD .and. njet .ne. nQCD-1) then |
1119 | - passcuts_user=.false. |
1120 | - return |
1121 | - endif |
1122 | - endif |
1123 | - 122 continue |
1124 | -c |
1125 | -c PHOTON (ISOLATION) CUTS |
1126 | -c |
1127 | -c find the photons |
1128 | + passcuts_jets=.false. |
1129 | + return |
1130 | + endif |
1131 | + endif |
1132 | + |
1133 | + return |
1134 | + end |
1135 | + |
1136 | + |
1137 | + |
1138 | + logical function passcuts_unlops_jv(p,istatus,ipdg,pQCD,nQCD,ickkw) |
1139 | + implicit none |
1140 | + include 'nexternal.inc' |
1141 | + integer istatus(nexternal) |
1142 | + integer iPDG(nexternal) |
1143 | + double precision p(0:4,nexternal) |
1144 | + logical is_a_j(nexternal) |
1145 | + integer nQCD |
1146 | + double precision pQCD(0:3,nexternal) |
1147 | + integer ickkw |
1148 | + include "cuts.inc" |
1149 | + double precision p_unlops(0:3,nexternal) |
1150 | + logical passUNLOPScuts |
1151 | + integer i, j |
1152 | + |
1153 | + REAL*8 pt |
1154 | + external pt |
1155 | + |
1156 | + passcuts_unlops_jv=.true. |
1157 | + |
1158 | + |
1159 | +c THE UNLOPS CUT: |
1160 | + if (ickkw.eq.4 .and. ptj.gt.0d0) then |
1161 | +c Use special pythia pt cut for minimal pT |
1162 | + do i=1,nexternal |
1163 | + do j=0,3 |
1164 | + p_unlops(j,i)=p(j,i) |
1165 | + enddo |
1166 | + enddo |
1167 | + call pythia_UNLOPS(p_unlops,passUNLOPScuts) |
1168 | + if (.not. passUNLOPScuts) then |
1169 | + passcuts_unlops_jv=.false. |
1170 | + return |
1171 | + endif |
1172 | +c THE VETO XSEC CUT: |
1173 | + elseif (ickkw.eq.-1 .and. ptj.gt.0d0) then |
1174 | +c Use veto'ed Xsec for analytic NNLL resummation |
1175 | + if (nQCD.ne.1) then |
1176 | + write (*,*) 'ERROR: more than one QCD parton in '/ |
1177 | + $ /'this event in cuts.f. There should only be one' |
1178 | + stop |
1179 | + endif |
1180 | + if (pt(pQCD(0,1)) .gt. ptj) then |
1181 | + passcuts_unlops_jv=.false. |
1182 | + return |
1183 | + endif |
1184 | + endif |
1185 | + return |
1186 | + end |
1187 | + |
1188 | + |
1189 | + |
1190 | + logical function passcuts_leptons(p,istatus,ipdg,is_a_lp_reco,is_a_lm_reco) |
1191 | + implicit none |
1192 | + include 'nexternal.inc' |
1193 | + integer istatus(nexternal) |
1194 | + integer iPDG(nexternal) |
1195 | + double precision p(0:4,nexternal) |
1196 | + logical is_a_lp_reco(nexternal),is_a_lm_reco(nexternal) |
1197 | + |
1198 | + REAL*8 R2_04,invm2_04,pt_04,eta_04,pt,eta |
1199 | + external R2_04,invm2_04,pt_04,eta_04,pt,eta |
1200 | + integer i,j |
1201 | + |
1202 | + include 'cuts.inc' |
1203 | + |
1204 | + passcuts_leptons=.true. |
1205 | + |
1206 | +c |
1207 | +c CHARGED LEPTON CUTS |
1208 | +c |
1209 | +c find the charged leptons (also used in the photon isolation cuts below) |
1210 | do i=1,nexternal |
1211 | - if (istatus(i).eq.1 .and. ipdg(i).eq.22 .and. .not.gamma_is_j) then |
1212 | - is_a_ph(i)=.true. |
1213 | - else |
1214 | - is_a_ph(i)=.false. |
1215 | + if(istatus(i).eq.1 .and. |
1216 | + & (ipdg(i).eq.11 .or. ipdg(i).eq.13 .or. ipdg(i).eq.15)) then |
1217 | + is_a_lm_reco(i)=.true. |
1218 | + else |
1219 | + is_a_lm_reco(i)=.false. |
1220 | + endif |
1221 | + if(istatus(i).eq.1 .and. |
1222 | + & (ipdg(i).eq.-11 .or. ipdg(i).eq.-13 .or. ipdg(i).eq.-15)) then |
1223 | + is_a_lp_reco(i)=.true. |
1224 | + else |
1225 | + is_a_lp_reco(i)=.false. |
1226 | endif |
1227 | enddo |
1228 | - if (ptgmin.gt.0d0) then |
1229 | - nph=0 |
1230 | - do j=nincoming+1,nexternal |
1231 | - if (is_a_ph(j)) then |
1232 | - nph=nph+1 |
1233 | - do i=0,3 |
1234 | - pgamma(i,nph)=p(i,j) |
1235 | - enddo |
1236 | - endif |
1237 | - enddo |
1238 | - if(nph.eq.0)goto 444 |
1239 | - write(*,*) 'ERROR in cuts.f: photon isolation is not working' |
1240 | - $ // ' for mixed QED-QCD corrections' |
1241 | - stop 1 |
1242 | - |
1243 | - if(isoEM)then |
1244 | - nem=nph |
1245 | - do k=1,nem |
1246 | - do i=0,3 |
1247 | - pem(i,k)=pgamma(i,k) |
1248 | - enddo |
1249 | - enddo |
1250 | - do j=nincoming+1,nexternal |
1251 | - if (is_a_lp(j).or.is_a_lm(j)) then |
1252 | - nem=nem+1 |
1253 | - do i=0,3 |
1254 | - pem(i,nem)=p(i,j) |
1255 | - enddo |
1256 | - endif |
1257 | - enddo |
1258 | - endif |
1259 | - |
1260 | - alliso=.true. |
1261 | - |
1262 | - j=0 |
1263 | - do while(j.lt.nph.and.alliso) |
1264 | -c Loop over all photons |
1265 | - j=j+1 |
1266 | - |
1267 | - ptg=pt(pgamma(0,j)) |
1268 | - if(ptg.lt.ptgmin)then |
1269 | - passcuts_user=.false. |
1270 | - return |
1271 | - endif |
1272 | - if (etagamma.gt.0d0) then |
1273 | - if (abs(eta(pgamma(0,j))).gt.etagamma) then |
1274 | - passcuts_user=.false. |
1275 | - return |
1276 | - endif |
1277 | - endif |
1278 | - |
1279 | -c Isolate from hadronic energy |
1280 | - do i=1,nQCD |
1281 | - drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pQCD(0,i))) |
1282 | - enddo |
1283 | - call sortzv(drlist,isorted,nQCD,ismode,isway,izero) |
1284 | - Etsum(0)=0.d0 |
1285 | - nin=0 |
1286 | - do i=1,nQCD |
1287 | - if(dble(drlist(isorted(i))).le.R0gamma)then |
1288 | - nin=nin+1 |
1289 | - Etsum(nin)=Etsum(nin-1)+pt(pQCD(0,isorted(i))) |
1290 | - endif |
1291 | - enddo |
1292 | - do i=1,nin |
1293 | - alliso=alliso .and. |
1294 | - $ Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))), |
1295 | - $ R0gamma,xn,epsgamma,ptg) |
1296 | - enddo |
1297 | - |
1298 | -c Isolate from EM energy |
1299 | - if(isoEM.and.nem.gt.1)then |
1300 | - do i=1,nem |
1301 | - drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pem(0,i))) |
1302 | - enddo |
1303 | - call sortzv(drlist,isorted,nem,ismode,isway,izero) |
1304 | -c First of list must be the photon: check this, and drop it |
1305 | - if(isorted(1).ne.j.or.drlist(isorted(1)).gt.1.e-4)then |
1306 | - write(*,*)'Error #1 in photon isolation' |
1307 | - write(*,*)j,isorted(1),drlist(isorted(1)) |
1308 | - stop |
1309 | - endif |
1310 | - Etsum(0)=0.d0 |
1311 | - nin=0 |
1312 | - do i=2,nem |
1313 | - if(dble(drlist(isorted(i))).le.R0gamma)then |
1314 | - nin=nin+1 |
1315 | - Etsum(nin)=Etsum(nin-1)+pt(pem(0,isorted(i))) |
1316 | +c apply the charged lepton cuts |
1317 | + do i=nincoming+1,nexternal |
1318 | + if (is_a_lp_reco(i).or.is_a_lm_reco(i)) then |
1319 | +c transverse momentum |
1320 | + if (ptl.gt.0d0) then |
1321 | + if (pt_04(p(0,i)).lt.ptl) then |
1322 | + passcuts_leptons=.false. |
1323 | + return |
1324 | + endif |
1325 | + endif |
1326 | +c pseudo-rapidity |
1327 | + if (etal.gt.0d0) then |
1328 | + if (abs(eta_04(p(0,i))).gt.etal) then |
1329 | + passcuts_leptons=.false. |
1330 | + return |
1331 | + endif |
1332 | + endif |
1333 | +c DeltaR and invariant mass cuts |
1334 | + if (is_a_lp_reco(i)) then |
1335 | + do j=nincoming+1,nexternal |
1336 | + if (is_a_lm_reco(j)) then |
1337 | + if (drll.gt.0d0) then |
1338 | + if (R2_04(p(0,i),p(0,j)).lt.drll**2) then |
1339 | + passcuts_leptons=.false. |
1340 | + return |
1341 | + endif |
1342 | + endif |
1343 | + if (mll.gt.0d0) then |
1344 | + if (invm2_04(p(0,i),p(0,j),1d0).lt.mll**2) then |
1345 | + passcuts_leptons=.false. |
1346 | + return |
1347 | + endif |
1348 | + endif |
1349 | + if (ipdg(i).eq.-ipdg(j)) then |
1350 | + if (drll_sf.gt.0d0) then |
1351 | + if (R2_04(p(0,i),p(0,j)).lt.drll_sf**2) then |
1352 | + passcuts_leptons=.false. |
1353 | + return |
1354 | + endif |
1355 | + endif |
1356 | + if (mll_sf.gt.0d0) then |
1357 | + if (invm2_04(p(0,i),p(0,j),1d0).lt.mll_sf**2) |
1358 | + $ then |
1359 | + passcuts_leptons=.false. |
1360 | + return |
1361 | + endif |
1362 | + endif |
1363 | + endif |
1364 | endif |
1365 | enddo |
1366 | - do i=1,nin |
1367 | - alliso=alliso .and. |
1368 | - $ Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))), |
1369 | - $ R0gamma,xn,epsgamma,ptg) |
1370 | - enddo |
1371 | endif |
1372 | -c End of loop over photons |
1373 | - enddo |
1374 | - if(.not.alliso)then |
1375 | - passcuts_user=.false. |
1376 | - return |
1377 | endif |
1378 | - 444 continue |
1379 | -c End photon isolation |
1380 | - endif |
1381 | + enddo |
1382 | + |
1383 | + return |
1384 | + end |
1385 | + |
1386 | + logical function passcuts_pdgs(p,istatus,ipdg) |
1387 | + implicit none |
1388 | + include 'nexternal.inc' |
1389 | + integer istatus(nexternal) |
1390 | + integer iPDG(nexternal) |
1391 | + double precision p(0:4,nexternal) |
1392 | +c PDG specific cut |
1393 | + double precision etmin(nincoming+1:nexternal-1) |
1394 | + double precision etmax(nincoming+1:nexternal-1) |
1395 | + double precision mxxmin(nincoming+1:nexternal-1,nincoming+1:nexternal-1) |
1396 | + common /to_cuts/etmin,etmax,mxxmin |
1397 | + REAL*8 invm2_04,pt_04 |
1398 | + external invm2_04,pt_04 |
1399 | +c temporary variable for caching locally computation |
1400 | + double precision tmpvar |
1401 | + integer i,j |
1402 | + |
1403 | + passcuts_pdgs = .true. |
1404 | + |
1405 | |
1406 | C |
1407 | C PDG SPECIFIC CUTS (PT/M_IJ) |
1408 | @@ -397,49 +666,24 @@ |
1409 | if(etmin(i).gt.0d0 .or. etmax(i).gt.0d0)then |
1410 | tmpvar = pt_04(p(0,i)) |
1411 | if (tmpvar.lt.etmin(i)) then |
1412 | - passcuts_user=.false. |
1413 | + passcuts_pdgs=.false. |
1414 | return |
1415 | elseif (tmpvar.gt.etmax(i) .and. etmax(i).gt.0d0) then |
1416 | - passcuts_user=.false. |
1417 | + passcuts_pdgs=.false. |
1418 | return |
1419 | endif |
1420 | endif |
1421 | do j=i+1, nexternal-1 |
1422 | if (mxxmin(i,j).gt.0d0)then |
1423 | if (invm2_04(p(0,i),p(0,j),1d0).lt.mxxmin(i,j)**2)then |
1424 | - passcuts_user=.false. |
1425 | + passcuts_pdgs=.false. |
1426 | return |
1427 | endif |
1428 | endif |
1429 | - enddo |
1430 | + enddo |
1431 | enddo |
1432 | - |
1433 | - |
1434 | -C*************************************************************** |
1435 | -C*************************************************************** |
1436 | -C PUT HERE YOUR USER-DEFINED CUTS |
1437 | -C*************************************************************** |
1438 | -C*************************************************************** |
1439 | -C |
1440 | -c$$$C EXAMPLE: cut on top quark pT |
1441 | -c$$$C Note that PDG specific cut are more optimised than simple user cut |
1442 | -c$$$ do i=1,nexternal ! loop over all external particles |
1443 | -c$$$ if (istatus(i).eq.1 ! final state particle |
1444 | -c$$$ & .and. abs(ipdg(i)).eq.6) then ! top quark |
1445 | -c$$$C apply the pT cut (pT should be large than 200 GeV for the event to |
1446 | -c$$$C pass cuts) |
1447 | -c$$$ if ( p(1,i)**2+p(2,i)**2 .lt. 200d0**2 ) then |
1448 | -c$$$C momenta do not pass cuts. Set passcuts_user to false and return |
1449 | -c$$$ passcuts_user=.false. |
1450 | -c$$$ return |
1451 | -c$$$ endif |
1452 | -c$$$ endif |
1453 | -c$$$ enddo |
1454 | -c |
1455 | - return |
1456 | - end |
1457 | - |
1458 | - |
1459 | + return |
1460 | + end |
1461 | |
1462 | |
1463 | C*************************************************************** |
1464 | @@ -983,5 +1227,21 @@ |
1465 | return |
1466 | end |
1467 | |
1468 | + integer function get_n_tagged_photons() |
1469 | + implicit none |
1470 | + integer i |
1471 | + include "nexternal.inc" |
1472 | + logical particle_tag(nexternal) |
1473 | + common /c_particle_tag/particle_tag |
1474 | + get_n_tagged_photons = 0 |
1475 | + |
1476 | + do i = nincoming+1, nexternal |
1477 | + if (particle_tag(i)) |
1478 | + $ get_n_tagged_photons = get_n_tagged_photons+1 |
1479 | + enddo |
1480 | + |
1481 | + return |
1482 | + end |
1483 | + |
1484 | |
1485 | |
1486 | |
1487 | === modified file 'Template/NLO/SubProcesses/fks_singular.f' |
1488 | --- Template/NLO/SubProcesses/fks_singular.f 2021-09-30 16:30:15 +0000 |
1489 | +++ Template/NLO/SubProcesses/fks_singular.f 2021-10-15 11:42:16 +0000 |
1490 | @@ -1519,6 +1519,12 @@ |
1491 | double precision wgt_ME_born,wgt_ME_real |
1492 | common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real |
1493 | |
1494 | + integer ntagph |
1495 | + double precision resc |
1496 | + integer get_n_tagged_photons |
1497 | + double precision get_rescale_alpha_factor |
1498 | + external get_n_tagged_photons get_rescale_alpha_factor |
1499 | + |
1500 | if (wgt1.eq.0d0 .and. wgt2.eq.0d0 .and. wgt3.eq.0d0) return |
1501 | c Check for NaN's and INF's. Simply skip the contribution |
1502 | if (wgt1.ne.wgt1) return |
1503 | @@ -1582,9 +1588,21 @@ |
1504 | icontr=icontr+1 |
1505 | call weight_lines_allocated(nexternal,icontr,max_wgt,max_iproc) |
1506 | itype(icontr)=type |
1507 | - wgt(1,icontr)=wgt1 |
1508 | - wgt(2,icontr)=wgt2 |
1509 | - wgt(3,icontr)=wgt3 |
1510 | + |
1511 | +C here we rescale the contributions by the ratio of alpha's in different |
1512 | +C schemes; it is needed when there are tagged photons around |
1513 | + ntagph = get_n_tagged_photons() |
1514 | + if (ntagph.eq.0) then |
1515 | + wgt(1,icontr)=wgt1 |
1516 | + wgt(2,icontr)=wgt2 |
1517 | + wgt(3,icontr)=wgt3 |
1518 | + else if (ntagph.gt.0) then |
1519 | + resc = get_rescale_alpha_factor(ntagph, orders(qed_pos)) |
1520 | + wgt(1,icontr) = wgt1 * resc |
1521 | + wgt(2,icontr) = wgt2 * resc |
1522 | + wgt(3,icontr) = wgt3 * resc |
1523 | + endif |
1524 | + |
1525 | bjx(1,icontr)=xbk(1) |
1526 | bjx(2,icontr)=xbk(2) |
1527 | scales2(1,icontr)=QES2 |
1528 | @@ -5555,6 +5573,8 @@ |
1529 | integer fks_j_from_i(nexternal,0:nexternal) |
1530 | & ,particle_type(nexternal),pdg_type(nexternal) |
1531 | common /c_fks_inc/fks_j_from_i,particle_type,pdg_type |
1532 | + logical particle_tag(nexternal) |
1533 | + common /c_particle_tag/particle_tag |
1534 | double precision particle_charge(nexternal) |
1535 | common /c_charges/particle_charge |
1536 | include "run.inc" |
1537 | @@ -5787,7 +5807,7 @@ |
1538 | if (particle_charge(i).eq.0d0.and.pdg_type(i).ne.22) |
1539 | $ cycle |
1540 | C set charge factors |
1541 | - if (pdg_type(i).eq.22) then |
1542 | + if (pdg_type(i).eq.22.and..not.particle_tag(i)) then |
1543 | c_used = 0d0 |
1544 | gamma_used = gamma_ph |
1545 | gammap_used = gammap_ph |
1546 | @@ -6067,7 +6087,7 @@ |
1547 | if (particle_charge(i).eq.0d0.and.pdg_type(i).ne.22) |
1548 | $ cycle |
1549 | C set charge factors |
1550 | - if (pdg_type(i).eq.22) then |
1551 | + if (pdg_type(i).eq.22.and..not.particle_tag(i)) then |
1552 | c_used = 0d0 |
1553 | gamma_used = gamma_ph |
1554 | gammap_used = gammap_ph |
1555 | @@ -6396,6 +6416,8 @@ |
1556 | double precision particle_charge(nexternal), particle_charge_born(nexternal-1) |
1557 | common /c_charges/particle_charge |
1558 | common /c_charges_born/particle_charge_born |
1559 | + logical particle_tag(nexternal) |
1560 | + common /c_particle_tag/particle_tag |
1561 | include 'coupl.inc' |
1562 | include 'q_es.inc' |
1563 | double precision p(0:3,nexternal),xmu2,double,single |
1564 | @@ -6504,7 +6526,7 @@ |
1565 | if (pdg_type(i).ne.22) then |
1566 | contr2=contr2-particle_charge(i)**2 |
1567 | contr1=contr1-3d0/2d0*particle_charge(i)**2 |
1568 | - else |
1569 | + elseif (.not.particle_tag(i)) then |
1570 | contr1=contr1-gamma_ph |
1571 | endif |
1572 | else |
1573 | @@ -6695,6 +6717,8 @@ |
1574 | double precision particle_charge(nexternal), particle_charge_born(nexternal-1) |
1575 | common /c_charges/particle_charge |
1576 | common /c_charges_born/particle_charge_born |
1577 | + logical particle_tag(nexternal) |
1578 | + common /c_particle_tag/particle_tag |
1579 | double precision zero |
1580 | parameter (zero=0d0) |
1581 | |
1582 | @@ -6839,7 +6863,7 @@ |
1583 | do i=nincoming+1,nexternal |
1584 | if (pdg_type(i).eq.21) ngluons_FKS(nFKSprocess) |
1585 | $ =ngluons_FKS(nFKSprocess)+1 |
1586 | - if (pdg_type(i).eq.22) nphotons_FKS(nFKSprocess) |
1587 | + if (pdg_type(i).eq.22.and..not.particle_tag(i)) nphotons_FKS(nFKSprocess) |
1588 | $ =nphotons_FKS(nFKSprocess)+1 |
1589 | enddo |
1590 | |
1591 | |
1592 | === modified file 'Template/NLO/SubProcesses/makefile_fks_dir' |
1593 | --- Template/NLO/SubProcesses/makefile_fks_dir 2020-11-27 13:41:46 +0000 |
1594 | +++ Template/NLO/SubProcesses/makefile_fks_dir 2021-10-15 11:42:16 +0000 |
1595 | @@ -32,6 +32,7 @@ |
1596 | $(patsubst %.f,%.o,$(wildcard b_sf_???.f)) \ |
1597 | born.o sborn_sf.o extra_cnt_wrapper.o \ |
1598 | $(patsubst %.f,%.o,$(wildcard born_cnt_*.f)) \ |
1599 | + rescale_alpha_tagged.o \ |
1600 | fks_Sij.o $(fastjetfortran_madfks) fks_singular.o \ |
1601 | montecarlocounter.o reweight_xsec.o boostwdir2.o \ |
1602 | initcluster.o cluster.o splitorders_stuff.o \ |
1603 | |
1604 | === modified file 'Template/NLO/SubProcesses/recmom.f' |
1605 | --- Template/NLO/SubProcesses/recmom.f 2020-09-16 12:55:48 +0000 |
1606 | +++ Template/NLO/SubProcesses/recmom.f 2021-10-15 11:42:16 +0000 |
1607 | @@ -1,12 +1,17 @@ |
1608 | - subroutine recombine_momenta(R, etaph, reco_l, reco_q, p_in, pdg_in, p_out, pdg_out) |
1609 | + subroutine recombine_momenta(R, etaph, reco_l, reco_q, p_in, |
1610 | + $pdg_in, is_nextph_iso, p_out, pdg_out, is_nextph_iso_reco) |
1611 | implicit none |
1612 | ! recombine photons with the closest fermion if the distance is |
1613 | ! less than R and if the rapidity of photons is < etaph (etaph < 0 |
1614 | ! means no cut). Output a new set of momenta and pdgs corresponding |
1615 | ! to the recombined particles. If recombination occurs the photon |
1616 | ! disappears from the output particles |
1617 | + ! If isolated photons exist (is_nextph_iso(i)=True), they are not |
1618 | + ! used for the recombination |
1619 | ! arguments |
1620 | include 'nexternal.inc' |
1621 | + include 'run.inc' |
1622 | + include 'cuts.inc' |
1623 | double precision R, etaph, p_in(0:4,nexternal), p_out(0:4,nexternal) |
1624 | logical reco_l, reco_q |
1625 | integer pdg_in(nexternal), pdg_out(nexternal) |
1626 | @@ -17,13 +22,24 @@ |
1627 | integer n_ph, i_ph |
1628 | integer i,j |
1629 | integer ifreco |
1630 | - double precision dreco, dthis |
1631 | + double precision dreco, dthis, iso_getdrv40 |
1632 | integer skip |
1633 | logical is_light_charged_fermion |
1634 | - double precision R2, eta |
1635 | + logical iso_in(nexternal), iso_out(nexternal) |
1636 | + double precision R2_04, eta_04 |
1637 | + ! use from photon isolation |
1638 | + logical is_nextph_iso(nexternal) |
1639 | +c integer nphiso |
1640 | +c double precision pgamma_iso(0:3,nexternal) |
1641 | +c common / photon / is_nextph_iso,nphiso,pgamma_iso |
1642 | + ! Save for reco photon isolation |
1643 | + logical is_nextph_iso_reco(nexternal) |
1644 | +c common / photon_reco / is_nextph_iso_reco |
1645 | ! |
1646 | integer times_reco |
1647 | common/to_times_reco/ times_reco |
1648 | + REAL*8 pt,eta |
1649 | + external pt,eta |
1650 | ! reset everything |
1651 | do j=1,nexternal |
1652 | pdg_out(j)=0 |
1653 | @@ -46,28 +62,36 @@ |
1654 | nq = 0 |
1655 | endif |
1656 | |
1657 | - ! count the photons |
1658 | + ! count the photons that are not isolated |
1659 | n_ph=0 |
1660 | do i=nincoming+1, nexternal |
1661 | - if (pdg_in(i).eq.id_ph.and. |
1662 | - $ (abs(eta(p_in(0,i))).lt.etaph.or.etaph.lt.0d0)) then |
1663 | + if (pdg_in(i).eq.id_ph .and. pt(p_in(0,i)).ne.0d0 .and. |
1664 | + $ (abs(eta_04(p_in(0,i))).lt.etaph.or.etaph.lt.0d0) |
1665 | + $ .and. .not.is_nextph_iso(i) ) then |
1666 | n_ph=n_ph+1 |
1667 | i_ph=i |
1668 | endif |
1669 | enddo |
1670 | + |
1671 | +c Define iso_in as the iso input from photon isolation |
1672 | + do i=nincoming+1,nexternal |
1673 | + iso_in(i)=is_nextph_iso(i) |
1674 | + enddo |
1675 | + |
1676 | if (n_ph.eq.0 .or. (nl.eq.0 .and. nq.eq.0)) then |
1677 | ! do nothing |
1678 | do j=1,nexternal |
1679 | pdg_out(j)=pdg_in(j) |
1680 | + iso_out(j)=iso_in(j) |
1681 | do i=0,4 |
1682 | p_out(i,j)=p_in(i,j) |
1683 | enddo |
1684 | enddo |
1685 | - return |
1686 | elseif (n_ph.eq.1) then |
1687 | ! do nothing for initial states |
1688 | do j=1,nincoming |
1689 | pdg_out(j)=pdg_in(j) |
1690 | + iso_out(j)=iso_in(j) |
1691 | do i=0,4 |
1692 | p_out(i,j)=p_in(i,j) |
1693 | enddo |
1694 | @@ -78,7 +102,7 @@ |
1695 | if (i_ph.gt.0) then |
1696 | do i = nincoming+1, nexternal |
1697 | if (is_light_charged_fermion(pdg_in(i),nq,nl)) then |
1698 | - dthis=dsqrt(R2(p_in(0,i_ph),p_in(0,i))) |
1699 | + dthis=dsqrt(R2_04(p_in(0,i_ph),p_in(0,i))) |
1700 | if (dthis.le.dreco) then |
1701 | dreco=dthis |
1702 | ifreco=i |
1703 | @@ -90,6 +114,7 @@ |
1704 | ! do nothing also for final states |
1705 | do j=nincoming+1,nexternal |
1706 | pdg_out(j)=pdg_in(j) |
1707 | + iso_out(j)=iso_in(j) |
1708 | do i=0,4 |
1709 | p_out(i,j)=p_in(i,j) |
1710 | enddo |
1711 | @@ -100,11 +125,13 @@ |
1712 | do j=nincoming+1,nexternal |
1713 | if (j.ne.i_ph.and.j.ne.ifreco) then |
1714 | pdg_out(j-skip)=pdg_in(j) |
1715 | + iso_out(j-skip)=iso_in(j) |
1716 | do i=0,4 |
1717 | p_out(i,j-skip)=p_in(i,j) |
1718 | enddo |
1719 | elseif (j.eq.ifreco) then |
1720 | pdg_out(j-skip)=pdg_in(j) |
1721 | + iso_out(j-skip)=iso_in(j) |
1722 | do i=0,3 |
1723 | p_out(i,j-skip)=p_in(i,j)+p_in(i,i_ph) |
1724 | enddo |
1725 | @@ -119,6 +146,35 @@ |
1726 | stop 1 |
1727 | endif |
1728 | |
1729 | + do i=nincoming+1,nexternal |
1730 | + is_nextph_iso_reco(i)=iso_out(i) |
1731 | + enddo |
1732 | + |
1733 | + return |
1734 | + end |
1735 | + |
1736 | + |
1737 | + |
1738 | + |
1739 | + |
1740 | + subroutine recombine_momenta_notagph(R, etaph, reco_l, reco_q, p_in, pdg_in, p_out, pdg_out) |
1741 | + implicit none |
1742 | + ! wraps recombine_momenta, it provides a function interface |
1743 | + ! for the case where no tagged photon exist |
1744 | + ! (mostly for backward complatibility) |
1745 | + ! arguments |
1746 | + include 'nexternal.inc' |
1747 | + double precision R, etaph, p_in(0:4,nexternal), p_out(0:4,nexternal) |
1748 | + logical reco_l, reco_q |
1749 | + integer pdg_in(nexternal), pdg_out(nexternal) |
1750 | + |
1751 | + logical is_iso_photon_in(nexternal), is_iso_photon_out(nexternal) |
1752 | + |
1753 | + is_iso_photon_in(:) = .false. |
1754 | + |
1755 | + call recombine_momenta(R, etaph, reco_l, reco_q, p_in, |
1756 | + $pdg_in, is_iso_photon_in, p_out, pdg_out, is_iso_photon_out) |
1757 | + |
1758 | return |
1759 | end |
1760 | |
1761 | |
1762 | === modified file 'UpdateNotes.txt' |
1763 | --- UpdateNotes.txt 2021-09-08 08:59:52 +0000 |
1764 | +++ UpdateNotes.txt 2021-10-15 11:42:16 +0000 |
1765 | @@ -1,6 +1,7 @@ |
1766 | Update notes for MadGraph5_aMC@NLO (in reverse time order) |
1767 | |
1768 | -3.1.2 (??/??/??) |
1769 | +3.2.1 (XX/XX/21) |
1770 | + DP+HS+IT+MZ: EW corrections can be computed for tagged photons |
1771 | MZ+CS: Fixes relevant to PineAPPL: |
1772 | - minor fix on the gluon pdg code in pineappl_interface.cc |
1773 | - the integration channels are threated the same way with or without PineAPPL |
1774 | |
1775 | === modified file 'madgraph/core/diagram_generation.py' |
1776 | --- madgraph/core/diagram_generation.py 2021-07-21 19:09:13 +0000 |
1777 | +++ madgraph/core/diagram_generation.py 2021-10-15 11:42:16 +0000 |
1778 | @@ -30,6 +30,7 @@ |
1779 | |
1780 | import madgraph.core.base_objects as base_objects |
1781 | import madgraph.various.misc as misc |
1782 | +import madgraph.fks.fks_tag as fks_tag |
1783 | from madgraph import InvalidCmd, MadGraph5Error |
1784 | from six.moves import range |
1785 | from six.moves import zip |
1786 | @@ -1713,6 +1714,15 @@ |
1787 | if leg['state'] == True] |
1788 | polids = [tuple(leg['polarization']) for leg in process_definition['legs'] \ |
1789 | if leg['state'] == True] |
1790 | + |
1791 | + # keep track of the 'is_tagged' property of the legs if needed |
1792 | + try: |
1793 | + fstags = [leg['is_tagged'] for leg in process_definition['legs'] \ |
1794 | + if 'is_tagged' in leg.keys() and leg['state'] == True] |
1795 | + |
1796 | + except KeyError: |
1797 | + fstags = [] |
1798 | + |
1799 | # Generate all combinations for the initial state |
1800 | for prod in itertools.product(*isids): |
1801 | islegs = [\ |
1802 | @@ -1735,10 +1745,16 @@ |
1803 | red_fsidlist.add(tuple(tag)) |
1804 | # Generate leg list for process |
1805 | leg_list = [copy.copy(leg) for leg in islegs] |
1806 | - leg_list.extend([\ |
1807 | - base_objects.Leg({'id':id, 'state': True, 'polarization': fslegs[i]['polarization']}) \ |
1808 | - for i,id in enumerate(prod)]) |
1809 | |
1810 | + if not fstags: |
1811 | + leg_list.extend([\ |
1812 | + base_objects.Leg({'id':id, 'state': True, 'polarization': fsleg['polarization']}) \ |
1813 | + for id, fsleg in zip(prod, fslegs)]) |
1814 | + else: |
1815 | + leg_list.extend([\ |
1816 | + fks_tag.TagLeg({'id':id, 'state': True, 'polarization': fsleg['polarization'], 'is_tagged': tag}) \ |
1817 | + for id, fsleg, tag in zip(prod, fslegs, fstags)]) |
1818 | + |
1819 | legs = base_objects.LegList(leg_list) |
1820 | |
1821 | # Check for crossed processes |
1822 | |
1823 | === modified file 'madgraph/fks/fks_base.py' |
1824 | --- madgraph/fks/fks_base.py 2021-05-23 14:17:24 +0000 |
1825 | +++ madgraph/fks/fks_base.py 2021-10-15 11:42:16 +0000 |
1826 | @@ -416,6 +416,7 @@ |
1827 | legs = [(leg.get('id'), leg) for leg in leglist] |
1828 | self.pdgs = array.array('i',[s[0] for s in legs]) |
1829 | self.colors = [leg['color'] for leg in leglist] |
1830 | + self.particle_tags = [leg['is_tagged'] for leg in leglist] |
1831 | if not self.process['perturbation_couplings'] == ['QCD']: |
1832 | self.charges = [leg['charge'] for leg in leglist] |
1833 | else: |
1834 | @@ -529,6 +530,13 @@ |
1835 | leg in self.born_amp['process']['legs']] |
1836 | |
1837 | |
1838 | + def get_is_tagged(self): |
1839 | + """return the list of the 'is_tagged' keys |
1840 | + of each leg in born_amp""" |
1841 | + return [leg.get('is_tagged') for \ |
1842 | + leg in self.born_amp['process']['legs']] |
1843 | + |
1844 | + |
1845 | def get_leglist(self): |
1846 | """return the leg list |
1847 | for the born amp""" |
1848 | |
1849 | === modified file 'madgraph/fks/fks_common.py' |
1850 | --- madgraph/fks/fks_common.py 2021-03-16 15:13:03 +0000 |
1851 | +++ madgraph/fks/fks_common.py 2021-10-15 11:42:16 +0000 |
1852 | @@ -263,7 +263,6 @@ |
1853 | if dict == {}: |
1854 | dict = find_pert_particles_interactions(model, pert) |
1855 | splittings = [] |
1856 | -#check that the leg is a qcd leg |
1857 | |
1858 | if leg.get('id') in dict['pert_particles']: |
1859 | part = model.get('particle_dict')[leg.get('id')] |
1860 | @@ -284,6 +283,13 @@ |
1861 | nsoft += 1 |
1862 | if nsoft >= 1: |
1863 | for split in split_leg(leg, parts, model): |
1864 | + # check if the leg is tagged, that |
1865 | + # the same particles appear also in the two daughters |
1866 | + if 'is_tagged' in leg.keys() and leg['is_tagged'] and \ |
1867 | + leg['id'] != split[0]['id'] and \ |
1868 | + leg['id'] != split[1]['id']: |
1869 | + continue |
1870 | + |
1871 | # add the splitting, but check if there is |
1872 | # an initial-state lepton if the flag |
1873 | # include_init_leptons is False |
1874 | @@ -833,6 +839,7 @@ |
1875 | -'charge', which gives the charge of the leg |
1876 | -'massless', boolean, true if leg is massless |
1877 | -'spin' which gives the spin of leg |
1878 | + -'is_tagged', boolean, true if leg is tagged in the final state |
1879 | -'is_part', boolean, true if leg is an particle |
1880 | -'self_antipart', boolean, true if leg is an self-conjugated particle |
1881 | """ |
1882 | @@ -846,13 +853,14 @@ |
1883 | self['charge'] = 0. |
1884 | self['massless'] = True |
1885 | self['spin'] = 0 |
1886 | + self['is_tagged'] = False |
1887 | self['is_part'] = True |
1888 | self['self_antipart'] = False |
1889 | |
1890 | def get_sorted_keys(self): |
1891 | """Return particle property names as a nicely sorted list.""" |
1892 | keys = super(FKSLeg, self).get_sorted_keys() |
1893 | - keys += ['fks', 'color','charge', 'massless', 'spin','is_part','self_antipart'] |
1894 | + keys += ['fks', 'color','charge', 'massless', 'spin','is_tagged','is_part','self_antipart',] |
1895 | return keys |
1896 | |
1897 | |
1898 | @@ -868,7 +876,7 @@ |
1899 | six.reraise(self.PhysicsObjectError, "%s is not a valid leg %s flag" % \ |
1900 | str(value), name) |
1901 | |
1902 | - if name in ['massless','self_antipart','is_part']: |
1903 | + if name in ['massless','is_tagged','self_antipart','is_part']: |
1904 | if not isinstance(value, bool): |
1905 | six.reraise(self.PhysicsObjectError, "%s is not a valid boolean for leg flag %s" % \ |
1906 | str(value), name) |
1907 | @@ -879,3 +887,4 @@ |
1908 | return super(FKSLeg,self).filter(name, value) |
1909 | |
1910 | |
1911 | + |
1912 | |
1913 | === modified file 'madgraph/fks/fks_helas_objects.py' |
1914 | --- madgraph/fks/fks_helas_objects.py 2021-02-08 22:58:11 +0000 |
1915 | +++ madgraph/fks/fks_helas_objects.py 2021-10-15 11:42:16 +0000 |
1916 | @@ -901,6 +901,7 @@ |
1917 | contains: |
1918 | -- colors |
1919 | -- charges |
1920 | + -- particle_tags |
1921 | -- i/j/ij fks, ij refers to the born leglist |
1922 | -- ijglu |
1923 | -- need_color_links |
1924 | @@ -919,6 +920,7 @@ |
1925 | if fksrealproc != None: |
1926 | self.isfinite = False |
1927 | self.colors = fksrealproc.colors |
1928 | + self.particle_tags = fksrealproc.particle_tags |
1929 | self.charges = fksrealproc.charges |
1930 | self.fks_infos = fksrealproc.fks_infos |
1931 | self.is_to_integrate = fksrealproc.is_to_integrate |
1932 | |
1933 | === added file 'madgraph/fks/fks_tag.py' |
1934 | --- madgraph/fks/fks_tag.py 1970-01-01 00:00:00 +0000 |
1935 | +++ madgraph/fks/fks_tag.py 2021-10-15 11:42:16 +0000 |
1936 | @@ -0,0 +1,75 @@ |
1937 | +################################################################################ |
1938 | +# |
1939 | +# Copyright (c) 2009 The MadGraph5_aMC@NLO Development team and Contributors |
1940 | +# |
1941 | +# This file is a part of the MadGraph5_aMC@NLO project, an application which |
1942 | +# automatically generates Feynman diagrams and matrix elements for arbitrary |
1943 | +# high-energy processes in the Standard Model and beyond. |
1944 | +# |
1945 | +# It is subject to the MadGraph5_aMC@NLO license which should accompany this |
1946 | +# distribution. |
1947 | +# |
1948 | +# For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch |
1949 | +# |
1950 | +################################################################################ |
1951 | + |
1952 | +"""Definitions of the objects needed both for MadFKS with tagged particles""" |
1953 | + |
1954 | +import madgraph.core.base_objects as MG |
1955 | + |
1956 | + |
1957 | +class MultiTagLeg(MG.MultiLeg): |
1958 | + """a daughter class of MultiLeg, with the extra possibility of specifying |
1959 | + whether a given leg is tagged or not, via the "is_tagged" key |
1960 | + """ |
1961 | + |
1962 | + def default_setup(self): |
1963 | + """Default values for all properties""" |
1964 | + super(MultiTagLeg, self).default_setup() |
1965 | + self['is_tagged'] = False |
1966 | + |
1967 | + def get_sorted_keys(self): |
1968 | + """Return particle property names as a nicely sorted list.""" |
1969 | + keys = super(MultiTagLeg, self).get_sorted_keys() |
1970 | + keys += ['is_tagged'] |
1971 | + return keys |
1972 | + |
1973 | + |
1974 | + def filter(self, name, value): |
1975 | + """Filter for valid leg property values.""" |
1976 | + |
1977 | + if name == 'is_tagged': |
1978 | + if not isinstance(value, bool): |
1979 | + raise self.PhysicsObjectError( \ |
1980 | + "%s is not a valid string for leg 'is_tagged' flag" \ |
1981 | + % str(value)) |
1982 | + return super(MultiTagLeg,self).filter(name, value) |
1983 | + |
1984 | + |
1985 | + |
1986 | +class TagLeg(MG.Leg): |
1987 | + """a daughter class of Leg, with the extra possibility of specifying |
1988 | + whether a given leg is tagged or not, via the "is_tagged" key |
1989 | + """ |
1990 | + |
1991 | + def default_setup(self): |
1992 | + """Default values for all properties""" |
1993 | + super(TagLeg, self).default_setup() |
1994 | + self['is_tagged'] = False |
1995 | + |
1996 | + def get_sorted_keys(self): |
1997 | + """Return particle property names as a nicely sorted list.""" |
1998 | + keys = super(TagLeg, self).get_sorted_keys() |
1999 | + keys += ['is_tagged'] |
2000 | + return keys |
2001 | + |
2002 | + |
2003 | + def filter(self, name, value): |
2004 | + """Filter for valid leg property values.""" |
2005 | + |
2006 | + if name == 'is_tagged': |
2007 | + if not isinstance(value, bool): |
2008 | + raise self.PhysicsObjectError( \ |
2009 | + "%s is not a valid string for leg 'is_tagged' flag" \ |
2010 | + % str(value)) |
2011 | + return super(TagLeg,self).filter(name, value) |
2012 | |
2013 | === modified file 'madgraph/interface/madgraph_interface.py' |
2014 | --- madgraph/interface/madgraph_interface.py 2021-08-20 21:40:18 +0000 |
2015 | +++ madgraph/interface/madgraph_interface.py 2021-10-15 11:42:16 +0000 |
2016 | @@ -93,6 +93,8 @@ |
2017 | import madgraph.various.misc as misc |
2018 | import madgraph.various.cluster as cluster |
2019 | |
2020 | +import madgraph.fks.fks_tag as fks_tag |
2021 | + |
2022 | import models as ufomodels |
2023 | import models.import_ufo as import_ufo |
2024 | import models.write_param_card as param_writer |
2025 | @@ -4927,6 +4929,21 @@ |
2026 | state = True |
2027 | continue |
2028 | |
2029 | + # check if the particle is tagged (!PART!) |
2030 | + if part_name.startswith('!') and part_name.endswith('!'): |
2031 | + part_name = part_name[1:-1] |
2032 | + is_tagged = True |
2033 | + elif part_name.endswith('!') and part_name.count('!') == 2 and part_name[:part_name.find('!')].isdigit(): |
2034 | + part_name = part_name.replace('!','') |
2035 | + is_tagged = True |
2036 | +# misc.sprint(part_name) |
2037 | + else: |
2038 | + is_tagged = False |
2039 | + |
2040 | + # check that only final-state particles are tagged |
2041 | + if is_tagged and not state: |
2042 | + raise self.InvalidCmd("initial particles cannot be tagged") |
2043 | + |
2044 | mylegids = [] |
2045 | polarization = [] |
2046 | if '{' in part_name: |
2047 | @@ -5006,6 +5023,9 @@ |
2048 | |
2049 | duplicate =1 |
2050 | if part_name in self._multiparticles: |
2051 | + # multiparticles cannot be tagged |
2052 | + if is_tagged: |
2053 | + raise self.InvalidCmd("Multiparticles cannot be tagged") |
2054 | if isinstance(self._multiparticles[part_name][0], list): |
2055 | raise self.InvalidCmd("Multiparticle %s is or-multiparticle" % part_name + \ |
2056 | " which can be used only for required s-channels") |
2057 | @@ -5036,12 +5056,26 @@ |
2058 | |
2059 | if mylegids: |
2060 | for _ in range(duplicate): |
2061 | - myleglist.append(base_objects.MultiLeg({'ids':mylegids, |
2062 | - 'state':state, |
2063 | - 'polarization': polarization})) |
2064 | + if LoopOption in ['virt','sqrvirt','tree','noborn']: |
2065 | + # check that no tagged particles exist in this mode |
2066 | + if is_tagged: |
2067 | + raise self.InvalidCmd( |
2068 | + "%s mode does not handle tagged particles" % LoopOption) |
2069 | + |
2070 | + myleglist.append(base_objects.MultiLeg({'ids':mylegids, |
2071 | + 'state':state, |
2072 | + 'polarization': polarization})) |
2073 | + else: |
2074 | + myleglist.append(fks_tag.MultiTagLeg({'ids':mylegids, |
2075 | + 'state':state, |
2076 | + 'polarization': polarization, |
2077 | + 'is_tagged':is_tagged})) |
2078 | else: |
2079 | raise self.InvalidCmd("No particle %s in model" % part_name) |
2080 | |
2081 | + if any(['is_tagged' in l.keys() and l['is_tagged'] for l in myleglist]): |
2082 | + logger.warning('The process involves tagged particles. Please consider citing arXiv:2106.02059 if relevant.') |
2083 | + |
2084 | # Apply the keyword 'all' for perturbed coupling orders. |
2085 | if perturbation_couplings.lower() in ['all', 'loonly']: |
2086 | if perturbation_couplings.lower() in ['loonly']: |
2087 | |
2088 | === modified file 'madgraph/iolibs/export_fks.py' |
2089 | --- madgraph/iolibs/export_fks.py 2021-07-21 19:12:36 +0000 |
2090 | +++ madgraph/iolibs/export_fks.py 2021-10-15 11:42:16 +0000 |
2091 | @@ -627,6 +627,16 @@ |
2092 | writers.FortranWriter(filename), |
2093 | matrix_element) |
2094 | |
2095 | + filename = 'a0Gmuconv.inc' |
2096 | + startfroma0 = self.write_a0gmuconv_file( |
2097 | + writers.FortranWriter(filename), |
2098 | + matrix_element) |
2099 | + |
2100 | + filename = 'rescale_alpha_tagged.f' |
2101 | + self.write_rescale_a0gmu_file( |
2102 | + writers.FortranWriter(filename), |
2103 | + startfroma0, matrix_element) |
2104 | + |
2105 | filename = 'orders.h' |
2106 | self.write_orders_c_header_file( |
2107 | writers.CPPWriter(filename), |
2108 | @@ -1096,6 +1106,23 @@ |
2109 | writer.writelines(text) |
2110 | |
2111 | |
2112 | + def write_a0gmuconv_file(self, writer, matrix_element): |
2113 | + """writes an include file with the informations about the |
2114 | + alpha0 < > gmu conversion, to be used when the process has |
2115 | + tagged photons |
2116 | + """ |
2117 | + |
2118 | + bool_dict = {True: '.true.', False: '.false.'} |
2119 | + bornproc = matrix_element.born_me['processes'][0] |
2120 | + startfromalpha0 = False |
2121 | + if any([l['is_tagged'] and l['id'] == 22 for l in bornproc['legs']]): |
2122 | + if 'loop_qcd_qed_sm_a0' in bornproc['model'].get('modelpath'): |
2123 | + startfromalpha0 = True |
2124 | + |
2125 | + text = 'logical startfroma0\nparameter (startfroma0=%s)\n' % bool_dict[startfromalpha0] |
2126 | + writer.writelines(text) |
2127 | + return startfromalpha0 |
2128 | + |
2129 | def write_orders_c_header_file(self, writer, amp_split_size, amp_split_size_born): |
2130 | """writes the header file including the amp_split_size declaration for amcblast |
2131 | """ |
2132 | @@ -1105,6 +1132,57 @@ |
2133 | writer.writelines(text) |
2134 | |
2135 | |
2136 | + def write_rescale_a0gmu_file(self, writer, startfroma0, matrix_element): |
2137 | + """writes the function that computes the rescaling factor needed in |
2138 | + the case of external photons |
2139 | + """ |
2140 | + |
2141 | + # get the model parameters |
2142 | + params = sum([v for v in self.model.get('parameters').values()], []) |
2143 | + parnames = [p.name.lower() for p in params] |
2144 | + |
2145 | + bornproc = matrix_element.born_me['processes'][0] |
2146 | + # this is to ensure compatibility with standard processes |
2147 | + if not any([l['is_tagged'] and l['id'] == 22 for l in bornproc['legs']]): |
2148 | + to_check = [] |
2149 | + expr = '1d0' |
2150 | + conv_pol = '0d0' |
2151 | + conv_fin = '0d0' |
2152 | + |
2153 | + elif startfroma0: |
2154 | + to_check = ['mdl_aewgmu', 'mdl_aew'] |
2155 | + base = 'mdl_aewgmu/mdl_aew' |
2156 | + exp = 'qed_pow/2d0-ntag' |
2157 | + expr = '(%s)**(%s)' % (base, exp) |
2158 | + conv_fin = '(qed_pow - ntagph * 2d0) * MDL_ECOUP_DGMUA0_UV_EW_FIN_ * born_wgt' |
2159 | + conv_pol = '(qed_pow - ntagph * 2d0) * MDL_ECOUP_DGMUA0_UV_EW_1EPS_ * born_wgt' |
2160 | + else: |
2161 | + to_check = ['mdl_aew', 'mdl_aew0'] |
2162 | + base = 'mdl_aew0/mdl_aew' |
2163 | + exp = 'ntag' |
2164 | + expr = '(%s)**(%s)' % (base, exp) |
2165 | + conv_fin = '- ntagph * 2d0 * MDL_ECOUP_DGMUA0_UV_EW_FIN_ * born_wgt' |
2166 | + conv_pol = '- ntagph * 2d0 * MDL_ECOUP_DGMUA0_UV_EW_1EPS_ * born_wgt' |
2167 | + |
2168 | + replace_dict = {'rescale_fact': expr, |
2169 | + 'virtual_a0Gmu_conv_finite': conv_fin, |
2170 | + 'virtual_a0Gmu_conv_pole': conv_pol} |
2171 | + |
2172 | + if not all(p in parnames for p in to_check): |
2173 | + raise fks_common.FKSProcessError( |
2174 | + 'Some parameters needed when there are tagged '+\ |
2175 | + 'photons cannot be found in the model.\n' +\ |
2176 | + 'Please load the correct model and restriction ' +\ |
2177 | + '(e.g loop_qcd_qed_sm_Gmu-a0 or loop_qcd_qed_sm_a0-Gmu)') |
2178 | + |
2179 | + file = open(os.path.join(_file_path, \ |
2180 | + 'iolibs/template_files/rescale_alpha_tagged.inc')).read() |
2181 | + file = file % replace_dict |
2182 | + |
2183 | + # Write the file |
2184 | + writer.writelines(file) |
2185 | + |
2186 | + |
2187 | |
2188 | def write_orders_file(self, writer, matrix_element): |
2189 | """writes the include file with the informations about coupling orders. |
2190 | @@ -2857,6 +2935,7 @@ |
2191 | col_lines = [] |
2192 | pdg_lines = [] |
2193 | charge_lines = [] |
2194 | + tag_lines = [] |
2195 | fks_j_from_i_lines = [] |
2196 | split_type_lines = [] |
2197 | for i, info in enumerate(fks_info_list): |
2198 | @@ -2870,6 +2949,9 @@ |
2199 | 'DATA (PARTICLE_CHARGE_D(%d, IPOS), IPOS=1, NEXTERNAL) / %s /'\ |
2200 | % (i + 1, ', '.join('%19.15fd0' % charg\ |
2201 | for charg in fksborn.real_processes[info['n_me']-1].charges) )) |
2202 | + tag_lines.append( \ |
2203 | + 'DATA (PARTICLE_TAG_D(%d, IPOS), IPOS=1, NEXTERNAL) / %s /' \ |
2204 | + % (i + 1, ', '.join(bool_dict[tag] for tag in fksborn.real_processes[info['n_me']-1].particle_tags) )) |
2205 | fks_j_from_i_lines.extend(self.get_fks_j_from_i_lines(fksborn.real_processes[info['n_me']-1],\ |
2206 | i + 1)) |
2207 | split_type_lines.append( \ |
2208 | @@ -2884,6 +2966,7 @@ |
2209 | pdgs = [l.get('id') for l in bornproc.get('legs')] + [-21] |
2210 | colors = [l.get('color') for l in bornproc.get('legs')] + [8] |
2211 | charges = [l.get('charge') for l in bornproc.get('legs')] + [0.] |
2212 | + tags = [l.get('is_tagged') for l in bornproc.get('legs')] + [False] |
2213 | |
2214 | fks_i = len(colors) |
2215 | # fist look for a colored legs (set j to 1 otherwise) |
2216 | @@ -2919,6 +3002,8 @@ |
2217 | % ', '.join([str(pdg) for pdg in pdgs])] |
2218 | charge_lines = ['DATA (PARTICLE_CHARGE_D(1, IPOS), IPOS=1, NEXTERNAL) / %s /' \ |
2219 | % ', '.join('%19.15fd0' % charg for charg in charges)] |
2220 | + tag_lines = ['DATA (PARTICLE_TAG_D(1, IPOS), IPOS=1, NEXTERNAL) / %s /' \ |
2221 | + % ', '.join(bool_dict[tag] for tag in tags)] |
2222 | fks_j_from_i_lines = ['DATA (FKS_J_FROM_I_D(1, %d, JPOS), JPOS = 0, 1) / 1, %d /' \ |
2223 | % (fks_i, fks_j)] |
2224 | split_type_lines = [ \ |
2225 | @@ -2930,6 +3015,7 @@ |
2226 | replace_dict['col_lines'] = '\n'.join(col_lines) |
2227 | replace_dict['pdg_lines'] = '\n'.join(pdg_lines) |
2228 | replace_dict['charge_lines'] = '\n'.join(charge_lines) |
2229 | + replace_dict['tag_lines'] = '\n'.join(tag_lines) |
2230 | replace_dict['fks_j_from_i_lines'] = '\n'.join(fks_j_from_i_lines) |
2231 | replace_dict['split_type_lines'] = '\n'.join(split_type_lines) |
2232 | |
2233 | |
2234 | === modified file 'madgraph/iolibs/template_files/fks_info.inc' |
2235 | --- madgraph/iolibs/template_files/fks_info.inc 2014-09-19 17:34:43 +0000 |
2236 | +++ madgraph/iolibs/template_files/fks_info.inc 2021-10-15 11:42:16 +0000 |
2237 | @@ -3,6 +3,7 @@ |
2238 | INTEGER EXTRA_CNT_D(%(nconfs)d), ISPLITORDER_BORN_D(%(nconfs)d), ISPLITORDER_CNT_D(%(nconfs)d) |
2239 | INTEGER FKS_J_FROM_I_D(%(nconfs)d, NEXTERNAL, 0:NEXTERNAL) |
2240 | INTEGER PARTICLE_TYPE_D(%(nconfs)d, NEXTERNAL), PDG_TYPE_D(%(nconfs)d, NEXTERNAL) |
2241 | + LOGICAL PARTICLE_TAG_D(%(nconfs)d, NEXTERNAL) |
2242 | REAL*8 PARTICLE_CHARGE_D(%(nconfs)d, NEXTERNAL) |
2243 | LOGICAL SPLIT_TYPE_D(%(nconfs)d, %(nsplitorders)d) |
2244 | LOGICAL NEED_COLOR_LINKS_D(%(nconfs)d), NEED_CHARGE_LINKS_D(%(nconfs)d) |
2245 | @@ -40,3 +41,9 @@ |
2246 | C Particle charge: |
2247 | C charge is set 0. with QCD corrections, which is irrelevant |
2248 | %(charge_lines)s |
2249 | + |
2250 | +C |
2251 | +C Tagged particles: |
2252 | +C |
2253 | +%(tag_lines)s |
2254 | + |
2255 | |
2256 | === added file 'madgraph/iolibs/template_files/rescale_alpha_tagged.inc' |
2257 | --- madgraph/iolibs/template_files/rescale_alpha_tagged.inc 1970-01-01 00:00:00 +0000 |
2258 | +++ madgraph/iolibs/template_files/rescale_alpha_tagged.inc 2021-10-15 11:42:16 +0000 |
2259 | @@ -0,0 +1,40 @@ |
2260 | +double precision function get_rescale_alpha_factor(ntag,qed_pow) |
2261 | +C returns the power of ratios of alpha needed to compensate |
2262 | +C for the presence of external tagged photons |
2263 | +C It is automatically written, and it detects whether the |
2264 | +C starting model is in the Gmu or alpha0 scheme. |
2265 | +C ntag is the number of tagged photons, qed_pow the power of gweak |
2266 | +C of the contribution considered |
2267 | +implicit none |
2268 | +integer ntag, qed_pow |
2269 | +INCLUDE '../../Source/MODEL/input.inc' |
2270 | + |
2271 | +get_rescale_alpha_factor = 1d0 |
2272 | +if (ntag.eq.0) return |
2273 | + |
2274 | +get_rescale_alpha_factor = %(rescale_fact)s |
2275 | + |
2276 | +return |
2277 | +end |
2278 | + |
2279 | + |
2280 | +double precision function get_virtual_a0Gmu_conv(qed_pow, ntagph, ivirt, born_wgt) |
2281 | +C returns the piece to compensate the a0<>Gmu change of scheme for the |
2282 | +C virtuals (single pole if ivirt = 1, finite if ivirt = 0 |
2283 | +implicit none |
2284 | +integer qed_pow, ntagph, ivirt |
2285 | +double precision born_wgt |
2286 | +INCLUDE '../../Source/MODEL/input.inc' |
2287 | + |
2288 | +if (ivirt.eq.1) then |
2289 | +! single |
2290 | +get_virtual_a0Gmu_conv = %(virtual_a0Gmu_conv_pole)s |
2291 | +else if (ivirt.eq.0) then |
2292 | +! finite part |
2293 | +get_virtual_a0Gmu_conv = %(virtual_a0Gmu_conv_finite)s |
2294 | +else |
2295 | +write(*,*) 'Error get_virtual_a0Gmu_conv: Invalid ivirt', ivirt |
2296 | +endif |
2297 | + |
2298 | +return |
2299 | +end |
2300 | |
2301 | === modified file 'madgraph/various/banner.py' |
2302 | --- madgraph/various/banner.py 2021-10-14 15:00:59 +0000 |
2303 | +++ madgraph/various/banner.py 2021-10-15 11:42:16 +0000 |
2304 | @@ -4528,6 +4528,7 @@ |
2305 | """Rules |
2306 | e+ e- beam -> lpp:0 ebeam:500 |
2307 | p p beam -> set maxjetflavor automatically |
2308 | + process with tagged photons -> gamma_is_j = false |
2309 | """ |
2310 | |
2311 | # check for beam_id |
2312 | @@ -4552,13 +4553,23 @@ |
2313 | if proc_characteristic['ninitial'] == 1: |
2314 | #remove all cut |
2315 | self.remove_all_cut() |
2316 | + |
2317 | + # check for tagged photons |
2318 | + tagged_particles = set() |
2319 | |
2320 | # Check if need matching |
2321 | min_particle = 99 |
2322 | max_particle = 0 |
2323 | for proc in proc_def: |
2324 | + for leg in proc['legs']: |
2325 | + if leg['is_tagged']: |
2326 | + tagged_particles.add(leg['id']) |
2327 | min_particle = min(len(proc['legs']), min_particle) |
2328 | max_particle = max(len(proc['legs']), max_particle) |
2329 | + |
2330 | + if 22 in tagged_particles: |
2331 | + self['gamma_is_j'] = False |
2332 | + |
2333 | matching = False |
2334 | if min_particle != max_particle: |
2335 | #take one of the process with min_particle |
2336 | |
2337 | === modified file 'tests/acceptance_tests/test_cmd_amcatnlo.py' |
2338 | --- tests/acceptance_tests/test_cmd_amcatnlo.py 2020-11-27 13:41:46 +0000 |
2339 | +++ tests/acceptance_tests/test_cmd_amcatnlo.py 2021-10-15 11:42:16 +0000 |
2340 | @@ -700,3 +700,38 @@ |
2341 | |
2342 | result = save_load_object.load_from_file('%s/HTML/results.pkl' % self.path) |
2343 | return result[run_name] |
2344 | + |
2345 | + |
2346 | + def test_generate_taggedph_nloew(self): |
2347 | + """test the param_card created is correct""" |
2348 | + |
2349 | + |
2350 | + text = """ |
2351 | + import model loop_qcd_qed_sm_Gmu-a0 |
2352 | + generate u u~ > !a! !a! [QED] |
2353 | + output %s |
2354 | + launch NLO |
2355 | + set lepphreco False |
2356 | + set quarkphreco False |
2357 | + """ % (self.path) |
2358 | + |
2359 | + interface = MGCmd.MasterCmd() |
2360 | + interface.no_notification() |
2361 | + |
2362 | + open(pjoin(self.tmpdir,'cmd'),'w').write(text) |
2363 | + |
2364 | + |
2365 | + os.system('rm -rf %s/RunWeb' % self.path) |
2366 | + os.system('rm -rf %s/Events/run_*' % self.path) |
2367 | + |
2368 | + interface.exec_cmd('import command %s' % pjoin(self.tmpdir, 'cmd')) |
2369 | + # test the lhe event file exists |
2370 | + self.assertTrue(os.path.exists('%s/Events/run_01/summary.txt' % self.path)) |
2371 | + self.assertTrue(os.path.exists('%s/Events/run_01/run_01_tag_1_banner.txt' % self.path)) |
2372 | + self.assertTrue(os.path.exists('%s/Events/run_01/res_0.txt' % self.path)) |
2373 | + self.assertTrue(os.path.exists('%s/Events/run_01/res_1.txt' % self.path)) |
2374 | + self.assertTrue(os.path.exists('%s/Events/run_01/alllogs_0.html' % self.path)) |
2375 | + self.assertTrue(os.path.exists('%s/Events/run_01/alllogs_1.html' % self.path)) |
2376 | + |
2377 | + check_html_page(self, pjoin(self.path, 'crossx.html')) |
2378 | + check_html_page(self, pjoin(self.path, 'HTML', 'run_01', 'results.html')) |
2379 | |
2380 | === added file 'tests/unit_tests/fks/test_fks_taggedphotons.py' |
2381 | --- tests/unit_tests/fks/test_fks_taggedphotons.py 1970-01-01 00:00:00 +0000 |
2382 | +++ tests/unit_tests/fks/test_fks_taggedphotons.py 2021-10-15 11:42:16 +0000 |
2383 | @@ -0,0 +1,80 @@ |
2384 | +############################################################################## |
2385 | +# |
2386 | +# Copyright (c) 2010 The MadGraph5_aMC@NLO Development team and Contributors |
2387 | +# |
2388 | +# This file is a part of the MadGraph5_aMC@NLO project, an application which |
2389 | +# automatically generates Feynman diagrams and matrix elements for arbitrary |
2390 | +# high-energy processes in the Standard Model and beyond. |
2391 | +# |
2392 | +# It is subject to the MadGraph5_aMC@NLO license which should accompany this |
2393 | +# distribution. |
2394 | +# |
2395 | +# For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch |
2396 | +# |
2397 | +################################################################################ |
2398 | +from __future__ import absolute_import |
2399 | +from __future__ import print_function |
2400 | +from cmd import Cmd |
2401 | +from six.moves import zip |
2402 | +""" Basic test of the command interface """ |
2403 | + |
2404 | +import unittest |
2405 | +import madgraph |
2406 | +import madgraph.interface.master_interface as mgcmd |
2407 | +import madgraph.interface.extended_cmd as ext_cmd |
2408 | +import madgraph.interface.amcatnlo_interface as amcatnlocmd |
2409 | +import os |
2410 | +import madgraph.fks.fks_helas_objects as fks_helas |
2411 | +import copy |
2412 | +import madgraph.iolibs.save_load_object as save_load_object |
2413 | +import madgraph |
2414 | +import madgraph.various.misc as misc |
2415 | +root_path = os.path.split(os.path.dirname(os.path.realpath( __file__ )))[0] |
2416 | +root_path = os.path.dirname(root_path) |
2417 | +# root_path is ./tests |
2418 | +pjoin = os.path.join |
2419 | + |
2420 | +class TestAMCatNLOEWTagPh(unittest.TestCase): |
2421 | + """ a suite of extra tests for the ew stuff """ |
2422 | + |
2423 | + def setUp(self): |
2424 | + self.interface = mgcmd.MasterCmd() |
2425 | + |
2426 | + def test_generate_fks_tagph(self): |
2427 | + """check that the generate command works as expected. |
2428 | + In particular the correct number of real-emission processes |
2429 | + and FKS partners""" |
2430 | + cmd_list = [ |
2431 | + 'u u~ > a a [real=QED]', |
2432 | + 'u u~ > !a! !a! [real=QED]', |
2433 | + 'u u~ > !2a! [real=QED]', |
2434 | + 'u u~ > 2!a! [real=QED]'] |
2435 | + |
2436 | + nrealproc_list = [9, 3, 3, 3] |
2437 | + fks_j_list = [[1,2,3], [1,2], [1,2], [1,2]] # 4 is not there in the first set for symmetry |
2438 | + |
2439 | + |
2440 | + for cmd, nrealproc, fks_j in \ |
2441 | + zip(cmd_list, nrealproc_list, fks_j_list): |
2442 | + self.interface.do_generate(cmd) |
2443 | + |
2444 | + fksprocess = self.interface._fks_multi_proc['born_processes'][0] |
2445 | + |
2446 | + self.assertEqual(len(fksprocess.real_amps), nrealproc) |
2447 | + self.assertEqual(set([r.fks_infos[0]['j'] for r in fksprocess.real_amps]), set(fks_j)) |
2448 | + |
2449 | + def test_invalid_syntax_tag(self): |
2450 | + |
2451 | + cmd = "u u~ > !a! !a!" |
2452 | + self.assertRaises(madgraph.InvalidCmd, self.interface.do_generate, cmd) |
2453 | + |
2454 | + cmd = "u u~ > !z! a [real=QED]" |
2455 | + self.assertRaises(madgraph.InvalidCmd, self.interface.do_generate, cmd) |
2456 | + |
2457 | + cmd = "u u~ > !t! t~ a [real=QED]" |
2458 | + self.assertRaises(madgraph.InvalidCmd, self.interface.do_generate, cmd) |
2459 | + |
2460 | + cmd = "u u~ > !a! a [real=QCD]" # is this valid ? I guess NOT [note this is QCD] |
2461 | + self.assertRaises(madgraph.InvalidCmd, self.interface.do_generate, cmd) |
2462 | + |
2463 | + |
Hi Marco,
I'll have a closer look next week. But for now:
- There is a text conflict in the UpdateNotes.txt
- Why do you add the new models? Shouldn't they part of FeynRules webpage so that they are simply downloaded when needed?
Best,
Rikkert