Heavy Higgs production with SM interferences

Dear Olivier,

Thanks a lot for such a prompt answer!

I tried to generate p p > e+ e-, and it worked like you described: I have negative weight events on my LHE.

For the process, I only have matrix1.f, which contains 4 diagrams. I've tried to return -|M|^2, and I have negative weight in my final LHE! That's great :)

I have reset everything to default, and only subtracted M2^2, M3^2 and M4^2, and I saw a negative weight event in the final LHE file, so I guess I'd done something wrong in my previous attempt!

I'll mark the answer as accepted since it seems to be OK. I'll reopen if something is wrong. In the meantime, I'll try to play with BW cutoff to see what happen :)

Thanks a lot for your help !

Cheers,
Sébastien

Thanks Olivier Mattelaer, that solved my question.

Hi,

  I would like to try out this model but I can't run it. I generated the UFO package with Mathematica 5 and feynrules 2.0 (latest). Then, when I try to generate events with MadGraph 5 v2.0 beta 3 I get two kind of errors:

MG5>generate g g > t t~ QS0=2 QCD=2 QED=2 / a z
Command "generate g g > t t~ QS0=2 QCD=2 QED=2 / a z" interrupted with error:
InvalidCmd : No particle qs0=2 in model

Then, if I change the input this way, I get another error:

MG5>generate g g > t t~ QCD=2 QED=2
Command "generate g g > t t~ QCD=2 QED=2" interrupted with error:
PhysicsObjectError : goldstone is not a valid property for this object: Particle

     Valid property are ['ghost', 'texname', 'is_part', 'name', 'self_antipart', 'color', 'LeptonNumber', 'width', 'charge', 'mass', 'counterterm', 'Y', 'antiname', 'line', 'propagating', 'spin', 'antitexname', 'pdg_code']
Please report this bug on https://bugs.launchpad.net/madgraph5
More information is found in 'MG5_debug'.
Please attach this file to your report.

It looks like a version mismatch. How do I solve it?

Cheers,
Riccardo

Oh, btw, I also tried the same UFO package with MadGraph v1.5.14. I get a different error:

mg5>generate g g > t t~ / a z QCD=2 QED=2 QS0=2
INFO: Trying process: g g > t t~ / a z QCD=2 QED=2 QS0=2

Command "generate g g > t t~ / a z QCD=2 QED=2 QS0=2 " interrupted in sub-command:
"generate g g > t t~ / a z QCD=2 QED=2 QS0=2" with error:
KeyError : '1'
Please report this bug on https://bugs.launchpad.net/madgraph5
More information is found in 'MG5_debug'.
Please attach this file to your report.

With some printout, the vertex is:

vertex = {
    'id': 87,
    'legs': [{
    'id': 21,
    'number': 1,
    'state': False,
    'from_group': True,
    'onshell': None
}, {
    'id': 21,
    'number': 2,
    'state': False,
    'from_group': True,
    'onshell': None
}, {
    'id': 6000049,
    'number': 1,
    'state': True,
    'from_group': True,
    'onshell': None
}]
}
couplings= {'1': 1}

order hierarchy = {'QED': 2, 'QO1': 8, 'QO0': 8, 'QS1': 8, 'QS0': 8, 'QCD': 1}

Does this mean anything to you experts? The "couplings" dictionary looks suspicious. For other vertices (eg id='68'), it look like
couplings= {'QED': 1, 'QO0': 1}

Thanks a lot,
R.

Hi,

for the first error,
did you correctly the letter/number after QS and/or is QS0 is a couplings type of your model (looks like it fails to recognize it as a coupling). Could you send me the model by email (<email address hidden>) so that I can check this point.

For the second error, this is FR/MadGraph compatibility issue.
This will be fixed in the release of next week.

For the third error, this is a typical problem when some coupling didn't have any coupling_order associate to those.
which are then associated by FR to couplings= {'1': 1}. Please contact the FR author to know how to fix this.

Cheers,

Olivier

Dear Olivier,

I am doing the same analysis now like Sebastien. This thread was a very good start for me, I had the same questions actually, and here it's very well explained. I would be curious now if this background subtraction can be done automatically. I have to run separately for different masses of A/H, different widths, for both A and H. Moreover, I specify the decays of the tops as well (because I learned from some discussions that when BSM particles decay to ttbar and we let Pythia do the decay of the tops, this might not be correctly done, with spins and BRs properly taken care of, like for SM ttbar). So in the end I have many subfolders of the type "P0_gg_" where I must change the matrix1.f. An automatic subtraction would be very useful.

This is my script.mg5 which generates background gg->tt + signal gg->H->tt (for mH=500GeV, widthH=10%=50GeV) and their interference. And I want to remove the background. I run this script with ./bin/mg5_aMC script.mg5 :

import model Higgs_Effective_Couplings_FormFact
define l+ = e+ mu+ ta+
define l- = e- mu- ta-

generate g g > t t~ / h1 QED=99 QCD=99, ( t > b w+, w+ > l+ vl), ( t~ > b~ w-, w- > j j)
add process g g > t t~ / h1 QED=99 QCD=99, ( t > b w+, w+ > j j), ( t~ > b~ w-, w- > l- vl~)
add process g g > t t~ / h1 QED=99 QCD=99, ( t > b w+, w+ > j j), ( t~ > b~ w-, w- > j j)
add process g g > t t~ / h1 QED=99 QCD=99, ( t > b w+, w+ > l+ vl), ( t~ > b~ w-, w- > l- vl~)

output ggtt500-w10p-H
launch -i
multi_run 400
set cut_decays F
set MT 172.5
set MH 500.6
set MP 10000
set WH 50.32
set YMTAU 0.68547
set YMT 349.5841
set YMB 1.6198
set nevents 100000
set ebeam1 4000
set ebeam2 4000

So my questions:

1) can I subtract the background automatically ?

2) or at least, to reduce the number of subfolders P0_gg[...] , could I skip mentioning the decays of the top and have just "generate g g > t t~ / h1 QED=99 QCD=99 " , and then let Pythia take care of the top decay when converting to EVNT and D3PDs and hoping for correct BRs ? Do I have to turn on MadSpin for this ?

Thank you,
Madalina

Dear MG authors,

Forget about my previous question, I managed it:

1) can I subtract the background automatically ?
-> it's no problem to do it manually .

2) or at least, to reduce the number of subfolders P0_gg[...] ,
-> I use a restrict file which diagonalizes the CKM matrix and sets m(c)=m(mu)=their Yukawa couplings=0. This reduces the no of P0_gg[...] and the total size to half

*****************************************************************

I would have another problem. For gg->A->ttbar, m(A)=500GeV, width=10%, I'm comparing m(ttbar) of signal+interference, on the tops immediately generated from gg fusion, no radiation involved, from the LHE file, between:

I) the DIRECT approach, when I subtract the background matrix element. Histogram is scaled with MG_xsec/Sum_of_weights

II) the INDIRECT approach, when I generate 2 types of files: a) signal+background+interference, b) background . I reconstruct m(ttbar) for both files , scale them with MG_xsec/Sum_of_weights , and subtract a-b .

Version: MadGraph5_aMC@NLO 5.2.0.1

I have to show you again the whole mg5 script, it gets again long, sorry about that, but it's necessary, thank you for your patience !!!

*****************************************************************

I) DIRECT APPROACH : 100k events

import model Higgs_Effective_Couplings_FormFact
define l+ = e+ mu+ ta+
define l- = e- mu- ta-

generate g g > t t~ / h QED=99 QCD=99, ( t > b w+, w+ > l+ vl), ( t~ > b~ w-, w- > j j)
add process g g > t t~ / h QED=99 QCD=99, ( t > b w+, w+ > j j), ( t~ > b~ w-, w- > l- vl~)
add process g g > t t~ / h QED=99 QCD=99, ( t > b w+, w+ > j j), ( t~ > b~ w-, w- > j j)
add process g g > t t~ / h QED=99 QCD=99, ( t > b w+, w+ > l+ vl), ( t~ > b~ w-, w- > l- vl~)

output DirectDir

Then I edit all matrix1,2,3,...,.f files:

After:

      MATRIX1 = 0.D0
      DO I = 1, NCOLOR
        ZTEMP = (0.D0,0.D0)
        DO J = 1, NCOLOR
          ZTEMP = ZTEMP + CF(J,I)*JAMP(J)
        ENDDO
        MATRIX1=MATRIX1+ZTEMP*DCONJG(JAMP(I))/DENOM(I)
      ENDDO

Comes:

      matrixbg=0d0
      DO I = 1, 2
        ZTEMP = (0.D0,0.D0)
        DO J = 1, 2
        ZTEMP = ZTEMP + CF(J,I)*JAMP(J)
        ENDDO
        MATRIXbg=MATRIXbg+ZTEMP*DCONJG(JAMP(I))/DENOM(I)
      ENDDO
      matrix1=matrix1-matrixbg

Then as usual:

      AMP2(1)=AMP2(1)+AMP(1)*DCONJG(AMP(1))
      AMP2(2)=AMP2(2)+AMP(2)*DCONJG(AMP(2))
      AMP2(3)=AMP2(3)+AMP(3)*DCONJG(AMP(3))
      AMP2(4)=AMP2(4)+AMP(4)*DCONJG(AMP(4))
      DO I = 1, NCOLOR
        JAMP2(I)=JAMP2(I)+JAMP(I)*DCONJG(JAMP(I))
      ENDDO

Then I continue my runs:

multi_run 1
set cut_decays F
set MT 172.5
set MH 10000
set MP 500
set WH1 49.63
set YMTAU 1.1351
set YMT 260.8092
set YMB 2.6825
set nevents 100000
set ebeam1 4000
set ebeam2 4000

*****************************************************************

II) INDIRECT APPROACH : same couplings, widths, masses as DIRECT approach.

a) SIGNAL+BACKGROUND+INTERFERENCE 40M events

import model Higgs_Effective_Couplings_FormFact
define l+ = e+ mu+ ta+
define l- = e- mu- ta-

generate g g > t t~ / h QED=99 QCD=99, ( t > b w+, w+ > l+ vl), ( t~ > b~ w-, w- > j j)
add process g g > t t~ / h QED=99 QCD=99, ( t > b w+, w+ > j j), ( t~ > b~ w-, w- > l- vl~)
add process g g > t t~ / h QED=99 QCD=99, ( t > b w+, w+ > j j), ( t~ > b~ w-, w- > j j)
add process g g > t t~ / h QED=99 QCD=99, ( t > b w+, w+ > l+ vl), ( t~ > b~ w-, w- > l- vl~)

output Indirect-SBI
launch -i
multi_run 400
set cut_decays F
set MT 172.5
set MH 10000
set MP 500
set WH1 49.63
set YMTAU 1.1351
set YMT 260.8092
set YMB 2.6825
set nevents 100000
set ebeam1 4000
set ebeam2 4000

b) SM BACKGROUND 40Mevents: same as II)a) but I remove "/ h QED=99 QCD=99"

Result : http://www.ifh.de/~stanescu/S-p-I-Dir-vs-Indir-amp2-non-zero/S-p-I-nsimpl-amp2-n0-decay-long-unw.png . NO MATCH. I use the unweighted LHE sample, no showering (only LHE level), no extra jets no restrict card. Result doesn't change much if I use the Weighted sample or the restrict card.

*****************************************************************

III) LET TOPS STABLE IN THE DIRECT APPROACH:

import model Higgs_Effective_Couplings_FormFact
define l+ = e+ mu+ ta+
define l- = e- mu- ta-

generate g g > t t~ / h QED=99 QCD=99

The rest stays the same. Indirect approach stays the same (so it also still specifies the decay of the tops). Now there is more or less MATCH: http://www.ifh.de/~stanescu/S-p-I-Dir-vs-Indir-amp2-non-zero/S-p-I-nsimpl-amp2-n0-decay-short-unw.png

*****************************************************************

CONCLUSION:

Would you maybe have an idea why:
- there is MATCH when Direct approach doesn't decay the tops, but Indirect approach decays the tops.
- there is NO MATCH when Direct approach decays the tops, and Indirect approach decays the tops.

Maybe I have to temper at combining the decay channels together, when I subtract the background matrix element.

I'd appreciate any small suggestion. My deepest regards,
Madalina

Another thing: to have the match, it would be tempting to let in MG only gg->ttbar, and let the tops be decayed later by Pythia, but then the spin correlations (therefore Phi angles between a fermion from top and equivalent from antitop) would not be maintained.

Bwt, I'm not using MadSpin.

Hi,

Concerning the fact that the direct method seems to give the expected result only when the decay is not present.
This probably means that you have an error with your implementation in presence of decay.

> Another thing: to have the match, it would be tempting to let in MG only
> gg->ttbar, and let the tops be decayed later by Pythia, but then the
> spin correlations (therefore Phi angles between a fermion from top and
> equivalent from anti top) would not be maintained.

Indeed that’s correct.

> Bwt, I'm not using MadSpin.

MadSpin is not going to help you except if you modify the matrix element used by madspin as well.

A third solution is to use the Bridge method, you do not have the full spin-correlation but you have the main effect of it.
This method is not part of MG5_aMC but I have a private version which is working (but without any nice user interface)

Cheers,

Olivier

On 26 Jan 2015, at 19:01, Madalina Stanescu-Bellu <email address hidden> wrote:

> Question #238699 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/238699
>
> Madalina Stanescu-Bellu posted a new comment:
> Another thing: to have the match, it would be tempting to let in MG only
> gg->ttbar, and let the tops be decayed later by Pythia, but then the
> spin correlations (therefore Phi angles between a fermion from top and
> equivalent from antitop) would not be maintained.
>
> Bwt, I'm not using MadSpin.
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Hi Olivier,

Thanks for your prompt answer, and for reading my long message !

>Concerning the fact that the direct method seems to give the expected result only when the decay is not present. This probably means that you have an error with your implementation in presence of decay.

Yes ! Would you have an idea what else should I do in p0_gg[...]/matrix.f , besides subtracting the background element like I explained above ? Or maybe something I should do outside P0_gg[...] channel directories, where the channels are combined ?

I'm also in discussions with Diogo Franzosi, he's the author of the model, we don't have any more clue how to proceed . And I don't know if this behavior depends on the model or on MG.

Cheers,
Madalina

Hi,

Actually, the only problem that I can see with your method,
is that you did not split the color flow of the photon and the 1/Nc part of the gluon propagator.

But such difference should have the same impact with and without decay.
Actually did you compare the case without decay with the built-in method of MG?
(we do not support interference with decay chain, at least not for the moment)

For the rest, I would check that indeed all matrix.f need the exact same subtraction.
(is it always the flow 1 and 2 that need to be remove/…).

Cheers,

Olivier

On 26 Jan 2015, at 21:11, Madalina Stanescu-Bellu <email address hidden> wrote:

> Question #238699 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/238699
>
> Madalina Stanescu-Bellu posted a new comment:
> Hi Olivier,
>
> Thanks for your prompt answer, and for reading my long message !
>
>> Concerning the fact that the direct method seems to give the expected
> result only when the decay is not present. This probably means that you
> have an error with your implementation in presence of decay.
>
> Yes ! Would you have an idea what else should I do in
> p0_gg[...]/matrix.f , besides subtracting the background element like I
> explained above ? Or maybe something I should do outside P0_gg[...]
> channel directories, where the channels are combined ?
>
> I'm also in discussions with Diogo Franzosi, he's the author of the
> model, we don't have any more clue how to proceed . And I don't know if
> this behavior depends on the model or on MG.
>
> Cheers,
> Madalina
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Dear Olivier,

I'm getting now a strange behaviour at the output number of events, and BRs. As a reminder, I run 100K events of gg->A->ttbar , m(A)=500GeV, where A=h1 and H=h, like this (our model has only the heavy neutral pseudoscala A and scalar H):

import model Higgs_Effective_Couplings_FormFact
define l+ = e+ mu+ ta+
define l- = e- mu- ta-

generate g g > t t~ / h QED=99 QCD=99, ( t > b w+, w+ > l+ vl), ( t~ > b~ w-, w- > j j)
add process g g > t t~ / h QED=99 QCD=99, ( t > b w+, w+ > j j), ( t~ > b~ w-, w- > l- vl~)
add process g g > t t~ / h QED=99 QCD=99, ( t > b w+, w+ > j j), ( t~ > b~ w-, w- > j j)
add process g g > t t~ / h QED=99 QCD=99, ( t > b w+, w+ > l+ vl), ( t~ > b~ w-, w- > l- vl~)

with removing the Matrix Element for the SM gg->tt. bwcutoff=15 . And I run that by setting AMP2 for the SM ttbar to 0, or by letting it unchanged (so !=0 ) .

I) AMP2!=0

the good news is the M(ttbar) and xsec of this Direct SignalPlusInterference approach matched very close to the Indirect case (with 30M of SM ttbar, 30M of SM ttbar + BSM gg->A->ttbar+ their interference, then subtraction of their root files): http://www.ifh.de/~stanescu/my-tmp/new/S-p-I-nsimpl-amp2-n0-decay-long-unw.png

the bad news is I get these bad BRs and improper no of output events :

Fully hadronic: 30%
Semileptonic: 45%
Dileptonic: 25%

56718 events, from which most are SM

II) AMP2==0

bad news is the M(ttbar) doesn't match: http://www.ifh.de/~stanescu/my-tmp/new/S-p-I-nsimpl-amp2-0-decay-long-unw.png

good news is the BRs are correct here:

Fully hadronic: 45%
Semileptonic: 45%
Dileptonic: 11%

100000 events, from which most are BSM

*********************************************************
QUESTIONS:
*********************************************************

1) would you know why I get bad BRs and output events no at I ) ? Something seems to create a divergence at the generation, so I get 56k events instead of 100k. Is it the bwcutoff ?

2) At least I understand why I) has mostly SM events, it's from what you explained above ( even though my bwcutoff is not infinite, but 15) : "On the other hand, if you keep AMP2(M2)!=2 then with infinite bwcutoff the ratio of events with S0 events written in file will be given by \int |M1|^2 / (|M1|^2 + |M2|^2) "

Thanks a lot,
Madalina

Hi,

Phase-space integration and generation of events are actually very complicated to perform.
The phase-space parametrization that you use is actually not optimized at all for such kind of function (i.e. interference).
This leads to some uneficiency and can explain why you fail to generate some events.

From your result it is quite clear to me that setting AMP2 to zero is a bad choice since it leads to part of the phase-space which are not probed at all.

Concerning the BR, I would guess that this might be linked to cut? Do you still have cut_decays =F ?

1) would you know why I get bad BRs and output events no at I ) ? Something seems to create a divergence at the generation, so I get 56k events instead of 100k. Is it the bwcutoff ?

Nothing to do with bwcutoff.

2) At least I understand why I) has mostly SM events, it's from what you explained above ( even though my bwcutoff is not infinite, but 15) : "On the other hand, if you keep AMP2(M2)!=2 then with infinite bwcutoff the ratio of events with S0 events written in file will be given by \int |M1|^2 / (|M1|^2 + |M2|^2) "

yes

Cheers,

Olivier

Dear Olivier,

I make a summary of what happens. So I generate signal+interference directly, for gg->A->ttbar, with tops decayed in MG (I specify all possible decays in the process card: fully hadronic, semileptonic, dileptonic), m(A)=500GeV, width=10%, by removing the SM gg->ttbar matrix element, with bwcutoff=15. I call this case "Direct". This I compare with the "Indirect " case> I run a few 10M events of background and a few 10M of signal+background+interference, then subtract their plots.

model: Higgs Effective Couplings Form Factor , with top loops , written by Diogo.

Results for "Direct" :

I. With AMP2=JAMP2=0 (for background):

a) m(ttbar) scaled with weights and xsec shows a smaller Xsec . But BR are ok (45% fully hadronic and semileptonic, 11% dileptonic), and all the kinematic distributions (m(ttbar), pT, eta, etc) look scalable to the "Indirect" distributions with constant values.

b) Weighted samples have weights only at 2 values : + val and -val, more at +val.

c) I asked for 100k events, I get 100k events in the unweighted sample.

Therefore, even though you said it shouldn't work, I'm tempted to choose this , because the other approach looks much worse:

II. I let AMP2 and JAMP2 as they are :

a) m(ttbar) matches the "Indirect" but BR are bad (20% dileptonic).

b) Weights in the weighted sample are again at +val1 and -val1, but the dileptonic events have an extra distribution after +val1 , therefore they get a BR 2x higher than normal (20%) : http://www.ifh.de/~stanescu/my-tmp/study-BR-weights-A/Direct-Amp2-n0/Direct-amp2-n0-Weighted-Weights-dilep.png . And this is fully hadronic: http://www.ifh.de/~stanescu/my-tmp/study-BR-weights-A/Direct-Amp2-n0/Direct-amp2-n0-Weighted-Weights-fully-hadr.png , with only 2 peaks.

c) Also, I asked for 100k events and I get 56k in the unweighted sample.

d) I managed to approximately fix the BR by setting artificially all weights to +1 (if they are +val) or -1 (for -val), and m(mu)=m(c)=0 , as well their Yukawa couplings .

Same behavior happens for the scalar H.

QUESTIONS:

1) Would you think it's worth trying in MadSpin ? You said I'd have to take care of removing the matrix element in Madspin too, I don't have time to figure that out and restart another study, because we have to hurry now to produce the samples, by end of next week.

2) where would MadSpin actually help me, to fix the fact that the narrow width approximation is not suitable for some kinds of decay patterns? But I noticed that when I generate signal gg->A/H->ttbar, width=1%, without MadSpin, I get good BR and other kinematics.

3) What if this increase BR for dileptonic to 20% is actually supposed to happen for interference ? But you said MG doesn't do interference at the decay products ?

4) I'll try quickly in HEFT model, the one without top loops , to see if I get same BR problem.

Thanks a lot,
Madalina

ANSWER

1) This is up to you, but I doubt that you will have something working in order to have result by the end of next week

2) I do not think so.

3) Not for the moment at least. Now that I think about it, it might be a good point since the NWA is a bad approximation for interference term.

4) ok

Cheers,

Olivier

ok, thanks a lot for your answers, especially during the weekend !!! I'm more clear now, but may ask please a bit more for clarification ?

a) so you're suggesting that, because MG does not perform interference on decay products, these BRs much different than normal for signal+interference come from a software/setup bug, instead of being a normal behavior for interference ?

b) and I don't understand yet where would MadSpin help/change my signal+interference pattern, can it look totally wrong because I didn't use MadSpin ?

Cheers & thanks,
Madalina

at b) by interference pattern I mean the shape of m(ttbar) and other kinematic distributions, I don't mean xsec or BR, that I can avoid, but I'm worried now that the shapes could be wrong because of not using MadSpin.

Hi Madalina,

> a) so you're suggesting that, because MG does not perform interference
> on decay products, these BRs much different than normal for
> signal+interference come from a software/setup bug, instead of being a
> normal behavior for interference ?

I do not know for sure, but this is a possibility.

> b) and I don't understand yet where would MadSpin help/change my
> signal+interference pattern, can it look totally wrong because I didn't
> use MadSpin ?

In principle MadSpin should always agree with MG since they basically rely on the same hypothesis.

Cheers,

Olivier

On 22 Mar 2015, at 19:31, Madalina Stanescu-Bellu <email address hidden> wrote:

> Question #238699 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/238699
>
> Madalina Stanescu-Bellu posted a new comment:
> ok, thanks a lot for your answers, especially during the weekend !!! I'm
> more clear now, but may ask please a bit more for clarification ?
>
> a) so you're suggesting that, because MG does not perform interference
> on decay products, these BRs much different than normal for
> signal+interference come from a software/setup bug, instead of being a
> normal behavior for interference ?
>
> b) and I don't understand yet where would MadSpin help/change my
> signal+interference pattern, can it look totally wrong because I didn't
> use MadSpin ?
>
> Cheers & thanks,
> Madalina
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Dear Olivier,

I think I finally fixed the problem. The samples Signal+Interference done in the "Direct" way (by subtracting the SM gg->ttbar contributions from the Matrix Element), where I let the AMP2 variable as it's calculated in matrix*.f , was always producing bad BRs (of ttbar and of W). But the kinematic distributions looked already matched to the "indirect" case, except 5% difference in the negative m(ttbar) range ("Indirect"=where I generated the 300M events of Background and 300M of Signal+Background+Interference).

On the other hand, when for "Direct" case I additionally set AMP2 of the background diagrams to 0, BRs were correct but all diagrams needed to be scaled (I wasn't sure if the scale for all possible kinematic was constant and the same, and also I didn't want to risk myself with this approach, because you said the xsec should not have changed, so I feared something was wrong).

I checked many things that could affects the cuts for the AMP2!=0 case, as you suggested, and nothing had much influence, but I finally realized that when I do one single run of 1M events (the maximum allowed size of 1 single run in MG) Signal+Interf, all ttbar BRs are perfect: fully hadronic 44.53%, semileptonic 44.49%, dileptonic 10.97%, (as opposed to 30%, 44%, 26% for 100k events) . BRs get more and more correct as I increased #Events. Also kinematics distributions still match.

Also , before, when I was asking for example for 200k events unweighted sample, MG was producing 60k (so 30%) , but now, when I ask 1M, is producing 750k (so 75%, the percentage improved).

This makes me realize it was all the time a problem of statistics. The "Direct"approach saves you from generating 300M events in the "Indirect" approach, but doing only 100k is still too little for convergence.

*********************************************

NOW THE PROBLEM IS : the 1M events sample is correct as a whole for the BRs, but if you take 100-200K subsamples of it, it will be again incorrect. And I can never know if the ATLAS production team or the analysis groups will take subsamples for their work instead of the whole. Plus for the ATLAS production team, a signal of 1M events might be still too big.

I've been suggested to do PREINTEGRATION. Would you have any idea about it ? To run a few millions events Signal+Interference, and then narrow it down to like 200k events.

Thanks a lot for your help so far , it was really useful !
My regards,
Madalina

Hi,

> NOW THE PROBLEM IS : the 1M events sample is correct as a whole for the
> BRs, but if you take 100-200K subsamples of it, it will be again
> incorrect.

Why is that? it should not be?
Are your events weighted?

> I've been suggested to do PREINTEGRATION. Would you have any idea about
> it ? To run a few millions events Signal+Interference, and then narrow
> it down to like 200k events.

Do you mean partial unweighting? (if your events are not un-weighted obviously)

Cheers,

Olivier

On 23 Apr 2015, at 13:41, Madalina Stanescu-Bellu <email address hidden> wrote:

> Question #238699 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/238699
>
> Madalina Stanescu-Bellu posted a new comment:
> Dear Olivier,
>
> I think I finally fixed the problem. The samples Signal+Interference
> done in the "Direct" way (by subtracting the SM gg->ttbar contributions
> from the Matrix Element), where I let the AMP2 variable as it's
> calculated in matrix*.f , was always producing bad BRs (of ttbar and of
> W). But the kinematic distributions looked already matched to the
> "indirect" case, except 5% difference in the negative m(ttbar) range
> ("Indirect"=where I generated the 300M events of Background and 300M of
> Signal+Background+Interference).
>
> On the other hand, when for "Direct" case I additionally set AMP2 of the
> background diagrams to 0, BRs were correct but all diagrams needed to be
> scaled (I wasn't sure if the scale for all possible kinematic was
> constant and the same, and also I didn't want to risk myself with this
> approach, because you said the xsec should not have changed, so I feared
> something was wrong).
>
> I checked many things that could affects the cuts for the AMP2!=0 case,
> as you suggested, and nothing had much influence, but I finally realized
> that when I do one single run of 1M events (the maximum allowed size of
> 1 single run in MG) Signal+Interf, all ttbar BRs are perfect: fully
> hadronic 44.53%, semileptonic 44.49%, dileptonic 10.97%, (as opposed to
> 30%, 44%, 26% for 100k events) . BRs get more and more correct as I
> increased #Events. Also kinematics distributions still match.
>
> Also , before, when I was asking for example for 200k events unweighted
> sample, MG was producing 60k (so 30%) , but now, when I ask 1M, is
> producing 750k (so 75%, the percentage improved).
>
> This makes me realize it was all the time a problem of statistics. The
> "Direct"approach saves you from generating 300M events in the "Indirect"
> approach, but doing only 100k is still too little for convergence.
>
> *********************************************
>
> NOW THE PROBLEM IS : the 1M events sample is correct as a whole for the
> BRs, but if you take 100-200K subsamples of it, it will be again
> incorrect. And I can never know if the ATLAS production team or the
> analysis groups will take subsamples for their work instead of the
> whole. Plus for the ATLAS production team, a signal of 1M events might
> be still too big.
>
> I've been suggested to do PREINTEGRATION. Would you have any idea about
> it ? To run a few millions events Signal+Interference, and then narrow
> it down to like 200k events.
>
>
> Thanks a lot for your help so far , it was really useful !
> My regards,
> Madalina
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Ah, yes, I took from that 1M events unweighted sample, the first 250k events, and the distributions look the same like the whole 1M sample. So you're right ! That's very good, so in worst case, I generate the 1M sample and then take only the first 250k from the LHE text file and throw away the others.

I've been told I could do it more elegant with this preintgration, which will only generate 250k. I asked the ATLAS MadGraph team and they sent me to read about Gridpack generation, I think this is it.

 I don't know what's partial unweighting,

Thanks again !!! Cheers.
Madalina

I meant running on the Gridpack and asking for a high enough accuracy, considering that the SM gg->ttbar that I subtracted from the MW, has a xsec 100x than the signal gg->A/H->ttbar which is left in the ME.

Hi

> I've been told I could do it more elegant with this preintgration, which
> will only generate 250k. I asked the ATLAS MadGraph team and they sent
> me to read about Gridpack generation, I think this is it.

Yes that’s a good idea.

Olivier

On 23 Apr 2015, at 15:51, Madalina Stanescu-Bellu <email address hidden> wrote:

> Question #238699 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/238699
>
> Madalina Stanescu-Bellu posted a new comment:
> Ah, yes, I took from that 1M events unweighted sample, the first 250k
> events, and the distributions look the same like the whole 1M sample. So
> you're right ! That's very good, so in worst case, I generate the 1M
> sample and then take only the first 250k from the LHE text file and
> throw away the others.
>
> I've been told I could do it more elegant with this preintgration, which
> will only generate 250k. I asked the ATLAS MadGraph team and they sent
> me to read about Gridpack generation, I think this is it.
>
> I don't know what's partial unweighting,
>
> Thanks again !!! Cheers.
> Madalina
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Dear MG authors,

I still have a minor question. My signal+interference production works, I have the samples and we're requesting the official ATLAS simulation now, and here I compare the m(ttbar) between

"Indirect" (black): with 300 Million generation separately of signal+background+interference, and 300 Million of background; then subtracting their root files

"Direct" (red/green) : 1 single run generation with the MG code changed to subtract the matrix element of the gg->ttbar background

for the positive interference the ratio is 1, but 1.05 for the negative interference. Would you know where this could come from?

http://www.ifh.de/~stanescu/my-tmp/Indirect-vs-Direct-tops-decayed-1M-official/1M-vs-250k/Indirect-vs-Direct-tops-decayed-unw-amp2-n0-1M-part-mlm-fastjet_020.M_ttbar_reco_from_parton_tops_II.png

I use the same setup, process definitions etc, when running "Indirect" and "Direct".

In the end I can consider the ratio 1, I don't need such a high precision. But if I'd have to introduce a k-factor , I'd have to generate those O(100) million events for each sample, in order to compare and get the scale factor.

Dear Madalina,

Do you use the default dynamical scale of MG? (a ckkw inspired scale)
If yes, this can explain such difference. Indeed the default scale is based on the Feynman diagram used for the generation of the events.
Therefore the scale used for your background sample and for your signal+background+interference background can be different for the same final state momentum.
On the other hand you expect the same between signal+background+interference and the direct method (depending of your method to do the subtraction actually)

In other word the indirect method would be slightly inconsistent but should be within scale uncertainty of the direct method.

If you use some simpler scale (easily available in 2.3.0) then I do not know.

Cheers,

Olivier

On 06 Aug 2015, at 11:37, Madalina Stanescu-Bellu <email address hidden> wrote:

> Question #238699 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/238699
>
> Madalina Stanescu-Bellu posted a new comment:
> Dear MG authors,
>
> I still have a minor question. My signal+interference production works,
> I have the samples and we're requesting the official ATLAS simulation
> now, and here I compare the m(ttbar) between
>
> "Indirect" (black): with 300 Million generation separately of
> signal+background+interference, and 300 Million of background; then
> subtracting their root files
>
> "Direct" (red/green) : 1 single run generation with the MG code changed
> to subtract the matrix element of the gg->ttbar background
>
> for the positive interference the ratio is 1, but 1.05 for the negative
> interference. Would you know where this could come from?
>
> http://www.ifh.de/~stanescu/my-tmp/Indirect-vs-Direct-tops-decayed-1M-
> official/1M-vs-250k/Indirect-vs-Direct-tops-decayed-unw-amp2-n0-1M-part-
> mlm-fastjet_020.M_ttbar_reco_from_parton_tops_II.png
>
> I use the same setup, process definitions etc, when running "Indirect"
> and "Direct".
>
> In the end I can consider the ratio 1, I don't need such a high
> precision. But if I'd have to introduce a k-factor , I'd have to
> generate those O(100) million events for each sample, in order to
> compare and get the scale factor.
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Hi Olivier,

I'm using MG5_aMC_v2_0_1 , and by default dynamical scale of MG do you mean pdfwgt ? I set it to false :

  .false. = pdfwgt ! for ickkw=1, perform pdf reweighting

ickkw I set it to 1 (so MLM), but actually in producing the official samples I set it to 0, we don't need jet matching / merging.

Here is the complete header of the Indirect case signal+background+interference, do you see anything there that could be the fault?

Thanks,
Madalina

http://www.ifh.de/~stanescu/my-tmp/header.txt

Hi,

I was speaking about those line:
  F = fixed_ren_scale ! if .true. use fixed ren scale
  F = fixed_fac_scale ! if .true. use fixed fac scale

and if you were using 2.3.0, you have one more interesting line:
-1 = dynamical_scale_choice ! Choose one of the preselected dynamical choices
with the different value defined here: https://answers.launchpad.net/mg5amcnlo/+faq/2014

So yes you are in the dynamical scale and therefore you use the default one that I describe in the previous email.
In that version of MG5, the only way to change the scale is to modify by hand the file SubProcesses/setscales.f

Cheers,

Olivier

On 06 Aug 2015, at 12:27, Madalina Stanescu-Bellu <email address hidden> wrote:

> Question #238699 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/238699
>
> Madalina Stanescu-Bellu posted a new comment:
> http://www.ifh.de/~stanescu/my-tmp/header.txt
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Hi again,

OK, I see now. So the scale factor is changed with the central M^2+pT^2, and from diagram to diagram ,so that's why the Direct and Indirect differ a bit (3%).

Does this happen only if ICKKW!=0 ? In those tests I accidentally let ICKKW=1 everywhere, even though I don't need it, because I don't have extra jets. My tops are decayed to W and b, and the W to j j or j lep , there is nowhere only j , to confuse the jet multiplicity. So in the official Signal+Interference production i set ICKKW=0.

But now I'd need a bit more instructions please: what should I change in my setup for comparing Direct vs Indirect, to make them match 100% ? Maybe this in the input script :

T = fixed_ren_scale ! if .true. use fixed ren scale
T = fixed_fac_scale ! if .true. use fixed fac scale

or something in setscale.f ? This is the file http://www.ifh.de/~stanescu/my-tmp/setscales.f , I see it returns 0 now but if I uncomment , it implements MT^2+pT^2.

Then the other question is for the official production of Direct S+I, it's correct how I do it now fixed_ren_scale=fixed_fac_scale=F, or should I set them to T, or change there setscale.f as well?

Thank you ,
Madalina

Hi,

> Does this happen only if ICKKW!=0 ? In those tests I accidentally let
> ICKKW=1 everywhere, even though I don't need it, because I don't have
> extra jets. My tops are decayed to W and b, and the W to j j or j lep ,
> there is nowhere only j , to confuse the jet multiplicity. So in the
> official Signal+Interference production i set ICKKW=0.

The procedure that I describe above is for the renormalisation scale both for ickkw=0 and ickkw=1,
the different for ickkw=1 is that each QCD vertex has it’s own scale, the QCD diagram consider for the scale choice depend of a kt-clustering
(or cone depending of your run_card).

> But now I'd need a bit more instructions please: what should I change in
> my setup for comparing Direct vs Indirect, to make them match 100% ?
> T = fixed_ren_scale ! if .true. use fixed ren scale
> T = fixed_fac_scale ! if .true. use fixed fac scale

Yes that is a way to make them match, but this is a set of parameter that is typically quite far from the
experimental curve since this is too crude.

> or something in setscale.f ? This is the file
> http://www.ifh.de/~stanescu/my-tmp/setscales.f , I see it returns 0 now
> but if I uncomment , it implements MT^2+pT^2.

Yes that’s much better.

> Then the other question is for the official production of Direct S+I,
> it's correct how I do it now fixed_ren_scale=fixed_fac_scale=F, or
> should I set them to T, or change there setscale.f as well?

That’s a question for your IT team.
If they use 2.3.0 for the official production, then you can choose the dynamical scheme from the run_card.

Cheers,

Olivier

On 06 Aug 2015, at 16:21, Madalina Stanescu-Bellu <email address hidden> wrote:

> Question #238699 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/238699
>
> Madalina Stanescu-Bellu posted a new comment:
> Hi again,
>
> OK, I see now. So the scale factor is changed with the central M^2+pT^2,
> and from diagram to diagram ,so that's why the Direct and Indirect
> differ a bit (3%).
>
> Does this happen only if ICKKW!=0 ? In those tests I accidentally let
> ICKKW=1 everywhere, even though I don't need it, because I don't have
> extra jets. My tops are decayed to W and b, and the W to j j or j lep ,
> there is nowhere only j , to confuse the jet multiplicity. So in the
> official Signal+Interference production i set ICKKW=0.
>
> But now I'd need a bit more instructions please: what should I change in
> my setup for comparing Direct vs Indirect, to make them match 100% ?
> Maybe this in the input script :
>
> T = fixed_ren_scale ! if .true. use fixed ren scale
> T = fixed_fac_scale ! if .true. use fixed fac scale
>
> or something in setscale.f ? This is the file
> http://www.ifh.de/~stanescu/my-tmp/setscales.f , I see it returns 0 now
> but if I uncomment , it implements MT^2+pT^2.
>
> Then the other question is for the official production of Direct S+I,
> it's correct how I do it now fixed_ren_scale=fixed_fac_scale=F, or
> should I set them to T, or change there setscale.f as well?
>
> Thank you ,
> Madalina
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Hi Olivier,

>> Then the other question is for the official production of Direct S+I,
>> it's correct how I do it now fixed_ren_scale=fixed_fac_scale=F, or
>> should I set them to T, or change there setscale.f as well?

> That’s a question for your IT team. If they use 2.3.0 for the official production, then you can choose the dynamical scheme from the run_card.

I'm using MG5_aMC_v2_0_1 . Is it wrong that I let the default values fixed_ren_scale=fixed_fac_scale=F, and setscales.f returns 0? I don't talk about matching with Indirect now, just if it's OK to leave the Direct like this . From your previous post I understood that the Direct is ok, just the Indirect a bit faulty.

Thanks,
Madalina

Hi,

In the indirect method you perform a subtraction with two different scales.
To my point of view, It is just weird to do that since it increase your theoretical uncertainty by two and forbid to have a nice agreement between the direct and indirect approach.
I would not use it due to that reason, but I do not think that this is really wrong.
So my suggestion is to modify by hand the setscales.f to use a pure kinematical scale (like HT or HT/2) which will reduce your theoretical uncertainty.

Cheers,

Olivier

On 06 Aug 2015, at 17:31, Madalina Stanescu-Bellu <email address hidden> wrote:

> Question #238699 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/238699
>
> Madalina Stanescu-Bellu posted a new comment:
> Hi Olivier,
>
>>> Then the other question is for the official production of Direct S+I,
>>> it's correct how I do it now fixed_ren_scale=fixed_fac_scale=F, or
>>> should I set them to T, or change there setscale.f as well?
>
>> That’s a question for your IT team. If they use 2.3.0 for the official
> production, then you can choose the dynamical scheme from the run_card.
>
> I'm using MG5_aMC_v2_0_1 . Is it wrong that I let the default values
> fixed_ren_scale=fixed_fac_scale=F, and setscales.f returns 0? I don't
> talk about matching with Indirect now, just if it's OK to leave the
> Direct like this . From your previous post I understood that the Direct
> is ok, just the Indirect a bit faulty.
>
> Thanks,
> Madalina
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Thanks a lot Oliver ! And I really appreciate the speed in your reply !

Cheers,
Madalina

Oliver -> Olivier

Hi Olivier,

Happy Easter ! We are back at the scales studies , so I would have one more question. When you have time , I would appreciate some help, but enjoy your holidays !

The questions: by HT you mean sum of pT , right? This I don't see yet implemented in the http://www.desy.de/~stanescu/my-tmp/setscales.f from my version 5.2.0.1, so I would do it like this :

     rscale=0d0
     do i=3,nexternal
      rscale=rscale+pt(P(0,i))
     enddo

     q2fact(1)=0d0
     do i=3,nexternal
      q2fact(1)= q2fact(1)+pt(P(0,i))
     enddo
     q2fact(2)=q2fact(1)

Looks OK ? By the way, my event has the tops decayed as well, so I think I should only add the pT of the tops, or only from its decay products.

And you recommended me to use HT instead of the total transverse energy of the event or sum of the transverse mass because HT is not dependend on the invariant mass?

Thank you,
Madalina

Dear Olivier,

Thank you for your previous explanations, in the meantime I modified setscales.f to return MT^2+pT^2 for both the factorization and renormalization scale, and I have, as before:

  F = fixed_ren_scale ! if .true. use fixed ren scale
  F = fixed_fac_scale ! if .true. use fixed fac scale
  91.1880 = scale ! fixed ren scale
  91.1880 = dsqrt_q2fact1 ! fixed fact scale for pdf1
  91.1880 = dsqrt_q2fact2 ! fixed fact scale for pdf2
  1 = scalefact ! scale factor for event-by-event scales

The match Direct-Indirect is indeed better.

Before , setscales.f was returning 0 for both the factorization and renormalization scale, and I had either ickkw=0 or ickkw=1, but I am a bit confused what dynamic scale did I have before, exactly, CKKM-based, as you said ? I ask just to be sure, since CKKM is a jet-matching and merging algorithm, and I do not shower inside MG, and these comparison Direct-Indirect are at parton level.

My regards,
Madalina

Ah, sorry, I forgot an sqrt for the renormalization scale , so my new renormalization scale is :
c--------------------------------------
c-- scale^2 = \sum_i (pt_i^2+m_i^2)
c--------------------------------------
      rscale=0d0
      do i=3,nexternal
       rscale=rscale+pt(P(0,i))**2+dot(p(0,i),p(0,i))
      enddo
      rscale=dsqrt(rscale)

and factorization scale :

c--------------------------------------
c-- scale^2 = \sum_i (pt_i^2+m_i^2)
c--------------------------------------
      q2fact(1)=0d0
      do i=3,nexternal
       q2fact(1)=q2fact(1)+pt(P(0,i))**2+dot(p(0,i),p(0,i))
      enddo
      q2fact(2)=q2fact(1)

But anyway, my question is only if the old dynamic scale was indeed CKKM.

Indeed CKKW is a matching/merging method, but in the matching part, a clustering of the jets are applied on the partonic events in order to perform an alpha_s and pdf reweighting. This is that procedure that we use in order to assign the scale.

Cheers,

Olivier

Hi Olivier, and how does MG chooses jets when I generate everything at parton level ? So no showring inside.

Cheers,
Madalina

Not sure to understand your question.
Do you mean which kind of clustering is it used? if yes this is a kt algorithm.

Cheers,

Olivier

I mean, I do not shower inside MG, everything is at parton level , the Pythia showering is not included in my generation, I take later the LHE and shower myself, or give the LHE to the ATLAS team to make full simulation.

So then how is a CKKM dynamic scale done with jets, in my MG simulation, if there are no jets ?

Thanks,
Madalina

We apply the clustering algorithm on the parton events.

Cheers,

Olivier

> On Apr 8, 2016, at 12:37, Madalina Stanescu-Bellu <email address hidden> wrote:
>
> Question #238699 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/238699
>
> Madalina Stanescu-Bellu posted a new comment:
> I mean, I do not shower inside MG, everything is at parton level , the
> Pythia showering is not included in my generation, I take later the LHE
> and shower myself, or give the LHE to the ATLAS team to make full
> simulation.
>
> So then how is a CKKM dynamic scale done with jets, in my MG simulation,
> if there are no jets ?
>
> Thanks,
> Madalina
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Hello,

I am just doing the same thing but for Run-2 (13 TeV) now. I have been working with Madalina closely but
I wanted to get your expert advice on the plot that I am sharing here. So I see a larger difference of about 10%-15% between
the direct and indirect approach for Run-2. I just normalize the histograms to the cross-section reported by MG. One thing I did notice that I am not using same pdfset as Run-1. I do specify in my script to use the 10800 pdf from lhadpf set but it ends up using nn23l0 instead. But I have made sure that all of my processes Direct as well as Indirect end up using same pdf set which in my case is the default one. Do you have any idea in mind as to why I could be seeing a larger difference than Run-1. ?

PS I am not sure how can I share/attach plot on this page?

Sincerely
Ruchika

Hi Ruchika,

>PS I am not sure how can I share/attach plot on this page?

It is not possible to attach file to questions. (Do not ask me why this is the platform of ubuntu).
The only way is to put it somewhere online where you can share it.

> One thing I did notice that I am not using same pdfset as Run-1. I do specify in my script to use the 10800 pdf from lhadpf set but it >ends up using nn23l0 instead.

You have to set the pdlabel to "lhapdf" if you want to use lhapdf set.

> I wanted to get your expert advice on the plot that I am sharing here.

I guess that you discuss those stuff first with Madalina and that I can consider that all the above point and other discussion on this topic do not cover this difference. Right?

Cheers,

Olivier

Dear MG authors,

I am doing the very same analysis like what Madalina did. We are using MG5 2.0.2 and the method you suggested here: use a pure kinematical scale by modifying the setscales.f by hand, to see if it can give us better agreement between the direct and indirect approach.

However, when we compared the Type-II 2HDM gg > H/A > ttbar + Interf results using default CKKM scale and pure kinematical scale, we found large cross section difference between them. For example, when using mass of pseudo scalar A:700GeV/c^2 and tanb:3.5, we got the M_ttbar distribution like this: https://www.dropbox.com/s/4grfkjnq4b6ou26/default_vs_dynamic.png?dl=0
You can see that the ratio is almost flat and the cross section for pure kinematical scale is about 111% as much as the cross section for default CKKM, which is not what we would expect.

It actually happened in all of our samples:
process xsec_CKKM xsec_{pT^2+M^2} xsec_{pT^2+M^2}/xsec_CKKM
ggHtt_SI_mH600wH128_tanb040 4.530e-01 5.023e-01 111%
ggAtt_SI_mA600wA180_tanb040 6.364e-01 7.405e-01 116%
ggHtt_SI_mH600wH082_tanb050 3.936e-01 4.372e-01 111%
ggAtt_SI_mA600wA115_tanb050 7.685e-01 8.823e-01 115%
ggHtt_SI_mH600wH025_tanb090 1.735e-01 1.924e-01 111%
ggAtt_SI_mA600wA035_tanb090 4.231e-01 4.800e-01 113%
ggHtt_SI_mH600wH009_tanb150 6.860e-02 7.604e-02 111%
ggAtt_SI_mA600wA012_tanb150 1.714e-01 1.939e-01 113%
ggHtt_SI_mH600wH001_tanb350 1.212e-02 1.348e-02 111%
ggAtt_SI_mA600wA002_tanb350 2.756e-02 3.125e-02 113%
ggHtt_SI_mH600wH000_tanb600 2.381e-03 2.722e-03 114%
ggAtt_SI_mA600wA001_tanb600 3.048e-03 3.703e-03 121%
ggHtt_SI_mH700wH170_tanb040 3.332e-01 3.692e-01 111%
ggAtt_SI_mA700wA214_tanb040 5.814e-01 6.710e-01 115%
ggHtt_SI_mH700wH109_tanb050 2.955e-01 3.273e-01 111%
ggAtt_SI_mA700wA137_tanb050 6.333e-01 7.216e-01 114%
ggHtt_SI_mH700wH055_tanb070 1.949e-01 2.155e-01 111%
ggAtt_SI_mA700wA070_tanb070 4.624e-01 5.232e-01 113%
ggHtt_SI_mH700wH009_tanb170 4.043e-02 4.463e-02 110%
ggAtt_SI_mA700wA011_tanb170 9.833e-02 1.109e-01 113%
ggHtt_SI_mH700wH003_tanb300 1.183e-02 1.310e-02 111%
ggAtt_SI_mA700wA003_tanb300 2.994e-02 3.374e-02 113%
ggHtt_SI_mH700wH033_tanb090 1.301e-01 1.437e-01 110%
ggAtt_SI_mA700wA042_tanb090 3.169e-01 3.577e-01 113%
ggHtt_SI_mH700wH012_tanb150 5.141e-02 5.665e-02 110%
ggAtt_SI_mA700wA015_tanb150 1.257e-01 1.418e-01 113%
ggHtt_SI_mH700wH002_tanb350 9.442e-03 1.042e-02 110%
ggAtt_SI_mA700wA002_tanb350 2.093e-02 2.364e-02 113%

The one using pure dynamical scale is always 11%~13% larger than the one using CKKM scale.

We suspect that it may be caused by the MG5 version differences, but we can't find available MG5 2.0.1 in any public website for further testing.

Would it be possible for you to give us this specific version for further testing?

Sincerely,
Yu-Heng

Hi,

> which is not what we would expect.

What do you expect?
A 15% difference by changing the scale seems very reasonable to me for LO computation.
Note that if you use the 2.5.0.beta version you will get directly the cross-section for the
> pT^2+M^2
> sqrts

> 1/2(pT^2+M^2)
> CKKW
and their associated error band (by multiplying and dividing by two)
(obtaining all those number takes an additional 20s on my laptop for 10k events generated)

For example, if I run the following in 2.5.0:
import model heft
generate g g > h > t t~
output GG_H_TT
launch -f

I have the following printout at the end of the run:

INFO: #***************************************************************************
#
# original cross-section: 0.00689593142743
# scale variation: +10.2% -8.87%
# central scheme variation: +12.6% -1.16e-09%
# PDF variation: +6.06% -6.06%
#
# dynamical scheme # 1 : 0.00728949 +9.98% - 9% # \sum ET
# dynamical scheme # 2 : 0.00704055 +10.3% -8.98% # \sum\sqrt{m^2+pt^2}
# dynamical scheme # 3 : 0.00776545 +9.77% -9.33% # 0.5 \sum\sqrt{m^2+pt^2}
# dynamical scheme # 4 : 0.00689593 +10.2% -8.87% # \sqrt{\hat s}
#***************************************************************************

So indeed for this process, the CKKW scheme seems to give the lowest cross-section compare to the other-scheme.
and difference of the order of 10% seems the typical variation related to different choice of scale.

> Would it be possible for you to give us this specific version for
> further testing?

you can create it yourself by running the following command
bzr branch lp:mg5amcnlo
cd mg5amcnlo
bzr revert -r 250

#obviously bzr revert -r 251 gives you 2.0.2 and so on.
# to have the relation between the version and the number you can look here:
http://bazaar.launchpad.net/~madteam/mg5amcnlo/series2.0/changes

Cheers,

Olivier

> On Sep 11, 2016, at 21:02, Yu-Heng Chen <email address hidden> wrote:
>
> Question #238699 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/238699
>
> Yu-Heng Chen posted a new comment:
> Dear MG authors,
>
> I am doing the very same analysis like what Madalina did. We are using
> MG5 2.0.2 and the method you suggested here: use a pure kinematical
> scale by modifying the setscales.f by hand, to see if it can give us
> better agreement between the direct and indirect approach.
>
> However, when we compared the Type-II 2HDM gg > H/A > ttbar + Interf results using default CKKM scale and pure kinematical scale, we found large cross section difference between them. For example, when using mass of pseudo scalar A:700GeV/c^2 and tanb:3.5, we got the M_ttbar distribution like this: https://www.dropbox.com/s/4grfkjnq4b6ou26/default_vs_dynamic.png?dl=0
> You can see that the ratio is almost flat and the cross section for pure kinematical scale is about 111% as much as the cross section for default CKKM, which is not what we would expect.
>
> It actually happened in all of our samples:
> process xsec_CKKM xsec_{pT^2+M^2} xsec_{pT^2+M^2}/xsec_CKKM
> ggHtt_SI_mH600wH128_tanb040 4.530e-01 5.023e-01 111%
> ggAtt_SI_mA600wA180_tanb040 6.364e-01 7.405e-01 116%
> ggHtt_SI_mH600wH082_tanb050 3.936e-01 4.372e-01 111%
> ggAtt_SI_mA600wA115_tanb050 7.685e-01 8.823e-01 115%
> ggHtt_SI_mH600wH025_tanb090 1.735e-01 1.924e-01 111%
> ggAtt_SI_mA600wA035_tanb090 4.231e-01 4.800e-01 113%
> ggHtt_SI_mH600wH009_tanb150 6.860e-02 7.604e-02 111%
> ggAtt_SI_mA600wA012_tanb150 1.714e-01 1.939e-01 113%
> ggHtt_SI_mH600wH001_tanb350 1.212e-02 1.348e-02 111%
> ggAtt_SI_mA600wA002_tanb350 2.756e-02 3.125e-02 113%
> ggHtt_SI_mH600wH000_tanb600 2.381e-03 2.722e-03 114%
> ggAtt_SI_mA600wA001_tanb600 3.048e-03 3.703e-03 121%
> ggHtt_SI_mH700wH170_tanb040 3.332e-01 3.692e-01 111%
> ggAtt_SI_mA700wA214_tanb040 5.814e-01 6.710e-01 115%
> ggHtt_SI_mH700wH109_tanb050 2.955e-01 3.273e-01 111%
> ggAtt_SI_mA700wA137_tanb050 6.333e-01 7.216e-01 114%
> ggHtt_SI_mH700wH055_tanb070 1.949e-01 2.155e-01 111%
> ggAtt_SI_mA700wA070_tanb070 4.624e-01 5.232e-01 113%
> ggHtt_SI_mH700wH009_tanb170 4.043e-02 4.463e-02 110%
> ggAtt_SI_mA700wA011_tanb170 9.833e-02 1.109e-01 113%
> ggHtt_SI_mH700wH003_tanb300 1.183e-02 1.310e-02 111%
> ggAtt_SI_mA700wA003_tanb300 2.994e-02 3.374e-02 113%
> ggHtt_SI_mH700wH033_tanb090 1.301e-01 1.437e-01 110%
> ggAtt_SI_mA700wA042_tanb090 3.169e-01 3.577e-01 113%
> ggHtt_SI_mH700wH012_tanb150 5.141e-02 5.665e-02 110%
> ggAtt_SI_mA700wA015_tanb150 1.257e-01 1.418e-01 113%
> ggHtt_SI_mH700wH002_tanb350 9.442e-03 1.042e-02 110%
> ggAtt_SI_mA700wA002_tanb350 2.093e-02 2.364e-02 113%
>
> The one using pure dynamical scale is always 11%~13% larger than the one
> using CKKM scale.
>
> We suspect that it may be caused by the MG5 version differences, but we
> can't find available MG5 2.0.1 in any public website for further
> testing.
>
> Would it be possible for you to give us this specific version for
> further testing?
>
> Sincerely,
> Yu-Heng
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Dear Olivier,

Thank you for your prompt reply and clear explanations.
It really saves us a lot of time not to do further check on this.

Sincerely,
Yu-Heng

Dear experts,

I am writing here as the topic discussed in this thread seems very relevant to what I am currently doing.

I am doing interference studies between vector-like quarks (Y->Wb) and the standard model (single top). This means that I generate the following processes with MadGraph:
generate p p > j w+ b~ / tp tp~ bp bp~ x x~ z h a QED=99 QCD=99
add process p p > j w- b / tp tp~ bp bp~ x x~ z h a QED=99 QCD=99

In addition I am removing the background in the MG simulations in the same way as described above: modifying matrix*.f files in order to remove the matrix element components corresponding to the background.

As a cross-check that I am not forgetting any background matrix element, I have generated the above processes, then modified the matrix*.f files to only generate background processes:

MATRIX1 = MATRIX_BKG

The resulting production I compared with a usual MG background production, i.e.:
generate p p > j w+ b~ / tp tp~ bp bp~ x x~ y y~ z h a QED=99 QCD=99
add process p p > j w- b / tp tp~ bp bp~ x x~ y y~ z h a QED=99 QCD=99

I found that the 2 productions match very well, when I set the AMP2 values related to the signal to 0:

https://www-zeuthen.desy.de/~dialexan/VLQ_SM_Interference/BkgComp_beforeCuts.eps

Before setting AMP2 to 0 I was observing some signal events in the production. This I find in contradiction to what is explained above, namely that the AMP2 values should not affect the result.

In addition I also find that when generating the background-subtracted interference pattern, the cure changes with respect to whether the AMP2 values of the background are set to 0 or not. See:

https://www-zeuthen.desy.de/~dialexan/VLQ_SM_Interference/IntYM900k024.pdf

I would start by asking if this could be clarified, i.e. what the AMP2 values are really doing and why including them or not is affecting my mass distributions. I think this could get me closer to understanding why my comparison between the DIRECT and INDIRECT (according to Madalina's definition above) distributions shows a mismatch:

https://www-zeuthen.desy.de/~dialexan/VLQ_SM_Interference/ IntMY700k095_comp.eps

(this plot was performed with AMP2!=0. No kinematic cuts were applied and bwcutoff=10000)

Many thanks in advance for your help.
Didier

Hi,

The Amp2 are just helping to define the weight of the multi-channelling.
They are just controlling how the code performs the multi-channel.
More details are available here:
hep-ph/0208156<http://arXiv.org/abs/hep-ph/0208156>

Cheers,

Olivier

On 21 Jul 2017, at 14:22, Didier Alexandre <<email address hidden><mailto:<email address hidden>>> wrote:

Question #238699 on MadGraph5_aMC@NLO changed:
https://answers.launchpad.net/mg5amcnlo/+question/238699

Didier Alexandre posted a new comment:
Dear experts,

I am writing here as the topic discussed in this thread seems very
relevant to what I am currently doing.

I am doing interference studies between vector-like quarks (Y->Wb) and the standard model (single top). This means that I generate the following processes with MadGraph:
generate p p > j w+ b~ / tp tp~ bp bp~ x x~ z h a QED=99 QCD=99
add process p p > j w- b / tp tp~ bp bp~ x x~ z h a QED=99 QCD=99

In addition I am removing the background in the MG simulations in the
same way as described above: modifying matrix*.f files in order to
remove the matrix element components corresponding to the background.

As a cross-check that I am not forgetting any background matrix element,
I have generated the above processes, then modified the matrix*.f files
to only generate background processes:

MATRIX1 = MATRIX_BKG

The resulting production I compared with a usual MG background production, i.e.:
generate p p > j w+ b~ / tp tp~ bp bp~ x x~ y y~ z h a QED=99 QCD=99
add process p p > j w- b / tp tp~ bp bp~ x x~ y y~ z h a QED=99 QCD=99

I found that the 2 productions match very well, when I set the AMP2
values related to the signal to 0:

https://www-
zeuthen.desy.de/~dialexan/VLQ_SM_Interference/BkgComp_beforeCuts.eps

Before setting AMP2 to 0 I was observing some signal events in the
production. This I find in contradiction to what is explained above,
namely that the AMP2 values should not affect the result.

In addition I also find that when generating the background-subtracted
interference pattern, the cure changes with respect to whether the AMP2
values of the background are set to 0 or not. See:

https://www-
zeuthen.desy.de/~dialexan/VLQ_SM_Interference/IntYM900k024.pdf

I would start by asking if this could be clarified, i.e. what the AMP2
values are really doing and why including them or not is affecting my
mass distributions. I think this could get me closer to understanding
why my comparison between the DIRECT and INDIRECT (according to
Madalina's definition above) distributions shows a mismatch:

https://www-zeuthen.desy.de/~dialexan/VLQ_SM_Interference/
IntMY700k095_comp.eps

(this plot was performed with AMP2!=0. No kinematic cuts were applied
and bwcutoff=10000)

Many thanks in advance for your help.
Didier

--
You received this question notification because you are an answer
contact for MadGraph5_aMC@NLO.

Dear Oliver,

Thanks a lot for your reply.
Ok so it is not clear to me at this stage why the setting of AMP2 change my invariant masses. Maybe the next step would be to verify if I am correctly subtracting the background matrix elements.

I would appreciate if you could take a look at the following snippet as an example:
https://www-zeuthen.desy.de/~dialexan/VLQ_SM_Interference/example.f

It contains the definition of the AMP related to the relevant diagrams. Then the matrix element (squared) is defined (MATRIX1) where the signal, background and interference terms are included. Then I include the definition of the matrix element terms related to the background only (MATRIX_BKG), and finally I do the subtraction : MATRIX1 = MATRIX1 - MATRIX_BKG:

Is there something I am missing in this procedure? One source of an inconsistency in the background definition I can think of would be the following:
Diagram 3, although a background diagram, contains a coupling terms from the VLQ model:

CALL FFV2_3_2(W(1,5),W(1,4),GC_42,GC_43,MDL_MY,MDL_WY,W(1,8))

where GC_42,GC_43 are couplings defined in the VLQ_UFO model

In a normal SM background generation, these terms would not be present. This is a wild guess, although I am not familiar at all with these functions and hence why a background diagram would contain a VLQ coupling term.

Many Thanks
Cheers,
Didier

Why not selecting the interference contribution in a similar way as
described in:

https://answers.launchpad.net/mg5amcnlo/+question/651188

which involves no fiddling with the fortran source code?

On Sat, Jul 22, 2017 at 10:23 AM, Didier Alexandre <
<email address hidden>> wrote:

> Question #238699 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/238699
>
> Didier Alexandre posted a new comment:
> Dear Oliver,
>
> Thanks a lot for your reply.
> Ok so it is not clear to me at this stage why the setting of AMP2 change
> my invariant masses. Maybe the next step would be to verify if I am
> correctly subtracting the background matrix elements.
>
> I would appreciate if you could take a look at the following snippet as an
> example:
> https://www-zeuthen.desy.de/~dialexan/VLQ_SM_Interference/example.f
>
> It contains the definition of the AMP related to the relevant diagrams.
> Then the matrix element (squared) is defined (MATRIX1) where the signal,
> background and interference terms are included. Then I include the
> definition of the matrix element terms related to the background only
> (MATRIX_BKG), and finally I do the subtraction : MATRIX1 = MATRIX1 -
> MATRIX_BKG:
>
> Is there something I am missing in this procedure? One source of an
> inconsistency in the background definition I can think of would be the
> following:
> Diagram 3, although a background diagram, contains a coupling terms from
> the VLQ model:
>
> CALL FFV2_3_2(W(1,5),W(1,4),GC_42,GC_43,MDL_MY,MDL_WY,W(1,8))
>
> where GC_42,GC_43 are couplings defined in the VLQ_UFO model
>
> In a normal SM background generation, these terms would not be present.
> This is a wild guess, although I am not familiar at all with these
> functions and hence why a background diagram would contain a VLQ
> coupling term.
>
> Many Thanks
> Cheers,
> Didier
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.
>

--
Valentin

Hi Valentin,

this looks much easier. I will try it. However, I noticed the model I am using (VLQ_UFO) does not have a coupling_order.py file, so this makes me wonder where madgraph reads that file from.

Thanks
Didier

Dear Didier,

Just a comment.
If you are going for the way suggested by Valentin, why not do it with FeynRules.
I think it is not only more legit both also more easy.

Sincerely,
Yu-Heng

Dear Valentin, Olivier

I used Velentin's method to generate VLQ & interference and it worked just fine. Interestingly, the results are in very good agreement to my previous method involving the matrix elements. See here:

https://www-zeuthen.desy.de/~dialexan/VLQ_SM_Interference/M700k095_int.eps

So, while this method with the new coupling order is very useful as it is much cleaner and less fiddly than the matrix element method, it is in good agreement with the latter, and hence in disagreement with the first method which consisted of subtracting the following two productions normalised to their respective cross- sections:

a) Signal + interference + SM background (single top + wjets)
b) SM background (single top + wjets)

For your reference, the comparison looks like this:
https://www-zeuthen.desy.de/~dialexan/VLQ_SM_Interference/IntMY700k095_comp.eps

Let me also point you to the two proc cards related to productions a) and b), as they are the only thing that changes (all other settings stay the same)

a) https://www-zeuthen.desy.de/~dialexan/VLQ_SM_Interference/proc_card_mg5_SIB.dat

b) https://www-zeuthen.desy.de/~dialexan/VLQ_SM_Interference/proc_card_mg5_Bkg.dat

I am out of ideas of what could be causing this discrepancy, so hopefully you can provide me with some clue.
Let me thank you at this stage for your very useful help so far.

Cheers
Didier

Hi,

How do you handle the scale used for alphas?
As mentioned above the default scale of madgraph is not appropriate for interference.
Especially when you make comparison since the different computation are using different running method.

Cheers,

Olivier
________________________________________
From: <email address hidden> <email address hidden> on behalf of Didier Alexandre <email address hidden>
Sent: Tuesday, July 25, 2017 7:33 PM
To: Olivier Mattelaer
Subject: Re: [Question #238699]: Heavy Higgs production with SM interferences

Question #238699 on MadGraph5_aMC@NLO changed:
https://answers.launchpad.net/mg5amcnlo/+question/238699

Didier Alexandre posted a new comment:
Dear Valentin, Olivier

I used Velentin's method to generate VLQ & interference and it worked
just fine. Interestingly, the results are in very good agreement to my
previous method involving the matrix elements. See here:

https://www-
zeuthen.desy.de/~dialexan/VLQ_SM_Interference/M700k095_int.eps

So, while this method with the new coupling order is very useful as it
is much cleaner and less fiddly than the matrix element method, it is in
good agreement with the latter, and hence in disagreement with the first
method which consisted of subtracting the following two productions
normalised to their respective cross- sections:

a) Signal + interference + SM background (single top + wjets)
b) SM background (single top + wjets)

For your reference, the comparison looks like this:
https://www-zeuthen.desy.de/~dialexan/VLQ_SM_Interference/IntMY700k095_comp.eps

Let me also point you to the two proc cards related to productions a)
and b), as they are the only thing that changes (all other settings stay
the same)

a) https://www-
zeuthen.desy.de/~dialexan/VLQ_SM_Interference/proc_card_mg5_SIB.dat

b) https://www-
zeuthen.desy.de/~dialexan/VLQ_SM_Interference/proc_card_mg5_Bkg.dat

I am out of ideas of what could be causing this discrepancy, so hopefully you can provide me with some clue.
Let me thank you at this stage for your very useful help so far.

Cheers
Didier

--
You received this question notification because you are an answer
contact for MadGraph5_aMC@NLO.

Hi Olivier,

ok, my dynamical scale parameter was set to default, now I changed it to option 2 (Total transverse energy of the event.) My comparison, although not perfect, improved:

https://www-zeuthen.desy.de/~dialexan/VLQ_SM_Interference/IntMY700k095_comp2.eps

I would say there is still a small slope pattern in the peak region (within +- 2%) while the tails show a larger mismatch.

Would you say this is within the scale uncertainty or could there still be another source of inconsistency between my productions? I am not sure how close I should expect my two methods to agree.

Thanks a lot
Didier

Please try with a fixed scale setting (i.e. set the option dynamical to
False in the run_card and specify its value there [leaving it to the
default Z mass is fine for the purpose of this comparison]).

With a fixed scale, you should find perfect agreement.

On Thu, Jul 27, 2017 at 12:03 AM, Didier Alexandre <
<email address hidden>> wrote:

> Question #238699 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/238699
>
> Didier Alexandre posted a new comment:
> Hi Olivier,
>
> ok, my dynamical scale parameter was set to default, now I changed it to
> option 2 (Total transverse energy of the event.) My comparison, although
> not perfect, improved:
>
> https://www-
> zeuthen.desy.de/~dialexan/VLQ_SM_Interference/IntMY700k095_comp2.eps
>
> I would say there is still a small slope pattern in the peak region
> (within +- 2%) while the tails show a larger mismatch.
>
> Would you say this is within the scale uncertainty or could there still
> be another source of inconsistency between my productions? I am not sure
> how close I should expect my two methods to agree.
>
> Thanks a lot
> Didier
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.
>

--
Valentin

Hi,

I will actually rather try with different seed and with larger statistical sample.
In any case, I would not worry too much about this since this is indeed much smaller (by huge factor) than scale uncertainty.
But this is obviously your call.

Cheers,

Olivier

Hi Valentin, Olivier,

@Valentin: ok, so I fixed both the renormalisation and factorisation to the Z mass, but still didn't get a better match between the curves, maybe suggesting there is still some additional inconsistency somewhere?

@Olivier: I don't think this is a problem of statistics as there is a slope at the peak (so it is a shape issue). For the subtraction I used in each case (S+I+B and B) 5 million events (500 jobs with a different random seed).

I think the current stage of my curves are sufficient for me to move on with the interference studies, but I would still be curious to understand if there are any additional details I am ignoring.

Thanks a lo
Didier