Jet matching/merging & DJR analysis with pythia8

Asked by Jack Y. Araz on 2016-07-26

Hi

I'm working on DJR plots from hepmc file generated by py8, to find the proper value for xqcut however I can not see 0 jet event domination separated from 1,2,3,4 jet events in DJR1 for instance. I'm having 0,1,2 jet events having more or less the same curve and 3,4 jet events are slightly separated. On the other hand when I do exactly the same thing with pythia6 I got exact same plots done in TASI Lectures 2013 (https://cp3.irmp.ucl.ac.be/projects/madgraph/attachment/wiki/TASISchool13/Tutorial_TASI2013.pdf) or done in arXiv:1401.6277 p. 193. Here is how I generate the events:

generate p p > l+ l- @0
add process p p > l+ l- j @1
add process p p > l+ l- j j @2
add process p p > l+ l- j j j @3
add process p p > l+ l- j j j j @4

(same with neutrinos but to get more precision I generated them separately and generated 25000 events each also tau included)

changes on run_card:
1 = ickkw
True = auto_ptj_mjj
50.0 = mmll
10.0 = xqcut

in pythia8 I used main89 with mlm comand file:

Main:numberOfEvents = 50000
Beams:frameType = 4

JetMatching:merge = on
JetMatching:scheme = 1
JetMatching:setMad = off (I also tried with on but nothing changed)
JetMatching:qCut = 10.0
JetMatching:coneRadius = 1.0
JetMatching:etaJetMax = 10.0
JetMatching:nJetMax = 4
JetMatching:doShowerKt = 0

Check:epTolErr = 1e-2

LHEFInputs:nSubruns = 2
Main:subrun = 0
Beams:LHEF = <l+l- events PATH>/unweighted_events.lhe
Main:subrun = 1
Beams:newLHEFsameInit = on
Main:LHEFskipInit = on
Beams:LHEF = <vl vl~ event PATH>/unweighted_events.lhe

in ma5:
./bin/ma5 -H
submset main.merging.check = true
submset main.merging.njets = 4
import <PATH>/Zlepneut.hepmc
submit

Note that I also tried to use 'set main.merging.ma5_mode = true' but it just shows sum and other's overlapping on the sum.

Also py8 generates only around ~14000 hepmc events where I'm giving 50000 event in total from lhe files.
Am I doing something wrong in pythia_card because there shouldn't be that much difference between py6 and py8 I believe.

Thank you very much
Best regards
Jack

Question information

Language:
English Edit question
Status:
Open
For:
MadGraph5_aMC@NLO Edit question
Assignee:
Stefan Prestel Edit question
Last query:
2016-08-11
Last reply:
2016-07-30

This question was reopened

Hi Jack,

This should be a problem of your PY8 analysis.
I therefore send this thread to PY8 expert.

Note that we are validating a MG5aMC-PY8 interface that should handle such merging automatically.
However this should not be release before at least a month.

Cheers,

Olivier

shibasipu (shibasipu) said : #2

Hi Olivier,

            >>Note that we are validating a MG5aMC-PY8 interface that should handle such merging automatically.
>>However this should not be release before at least a month.

it is indeed a great news for us. Like Pythia6(the one interfaced with Madgraph) does the showering and hadronization when we on Pythia, Is the same thing is also possible for MG5aMC-PY8 interface ?

Regards,
Shiba

It will be possible yes.

olivier
> On Jul 26, 2016, at 11:37, shibasipu <email address hidden> wrote:
>
> Question #303531 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/303531
>
> shibasipu posted a new comment:
> Hi Olivier,
>
>>> Note that we are validating a MG5aMC-PY8 interface that should handle such merging automatically.
>>> However this should not be release before at least a month.
>
> it is indeed a great news for us. Like Pythia6(the one interfaced with
> Madgraph) does the showering and hadronization when we on Pythia, Is the
> same thing is also possible for MG5aMC-PY8 interface ?
>
>
> Regards,
> Shiba
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

shibasipu (shibasipu) said : #4

Hi Olivier,
                    thanks a lot. Another thing, Will there be any limit on no. of events for MG5+PY8 for a single run like for MG5+PY6 maximum one can generate 50,000 events ?

Regards,
Shiba

Stefan Prestel (prestel) said : #5

Hi Jack,

As Olivier mentioned, it would be good to see your analysis to see what's going on. For Pythia6, the DJR plots you get from MG5 are produced internally, and the "tagging" of an event as an "njet-event" also relies on internal information - to have maximal discriminating power.

It's e.g. not quite clear to me how you would decide to call an event a "0 jet event". Maybe that would already explain the behaviour.

Regarding the number of events: Pythia8 interprets the 50K you request as the number of tried events (i.e. that can be read from LHEF). The 14K that you get out are the accepted events, after vetoing to remove double-counting has taken place. For Z+jets, an efficiency of ~30% can be expected. So I wouldn't worry about this.

Cheers,
Stefan

Jack Y. Araz (jackaraz) said : #6

Dear Stefan, Oliver

Thanks for the quick response and I'm eagerly waiting for that update thank you very much.

> It's e.g. not quite clear to me how you would decide to call an event a "0 jet event". Maybe that would already explain the behaviour.

If you are referring the name tagging @0,@1 etc I also tried it without them and the result was the same.

>As Olivier mentioned, it would be good to see your analysis to see what's going on.

I added lhe files with hepmc file to this link (https://drive.google.com/folderview?id=0B6ioFrUfVu51eWFPV3JDd0tPZzA&usp=sharing) note that lhe file are created without @0,@1 but rest is the same, hepmc file on the other hand is created with the name tagged lhe files since I'm constantly trying every bit of information that I could find, I didn't change it but if you need to see the resent ones I can add them too.

Thank you very much
Best regards
Jack

Jack Y. Araz (jackaraz) said : #7

Dear Stefan

I also tried to get DJR plots directly from py8 to see if the problem is with ma5 but I'm keep getting an error. I added this piece of script to main89:

#include "Pythia8Plugins/JetMatchingPy8Internal.h"
...
JetMatchingMadgraph *fJetMatchingPy8InternalHook = NULL;
fJetMatchingPy8InternalHook = new JetMatchingMadgraph;
pythia.setUserHooksPtr(fJetMatchingPy8InternalHook);
const std::vector<double> djrmatch = fJetMatchingPy8InternalHook->GetDJR();
std::ofstream myoutput;
myoutput.open("background/DJR.txt");
for(unsigned int i=0; i < std::min((float)djrmatch.size(), (float)4);i++)
    {
    myoutput << djrmatch[i]<< "\t"
             << fJetMatchingPy8InternalHook->nMEpartons().first <<"\t"
             << fJetMatchingPy8InternalHook->nMEpartons().second
             << std::endl;
    }

where JetMatchingPy8Internal.h is taken from this link (https://cmssdt.cern.ch/SDT/doxygen/CMSSW_7_1_3/doc/html/d1/d23/JetMatchingPy8Internal_8h_source.html#l00185) and during the compilation I'm keep getting following error:

main89_djr.cc:215:69: error: ‘class Pythia8::JetMatchingMadgraph’ has no member named ‘GetDJR’
   const std::vector<double> djrmatch = fJetMatchingPy8InternalHook->GetDJR();

where if my limited knowled on c++ didnt mislead me GetDJR is member in the class. Thus I thought this approach may solve it if there is an interface bug in ma5 for py8 hepmc output.

Thank you very much
Best regards
Jack

Stefan Prestel (prestel) said : #8

-Hi Jack,

First, sorry - I didn't realise that you were producing the DJRs with MadAnalysis. I unfortunately have no experience regarding that. Here is an option that you could try, minimally editing main89.cc. I'm not sure how useful that is for your issue, and maybe you should not waste too much time on this.

Now for your debuggong: It is probably for the best to remove the external code that you've linked. It's overcomplicating things, and I also cannot access the page. Instead, I suggest the following:

- In main89.cc, replace the line

UserHooks* matching = NULL;

with

JetMatchingMadgraph* matching = NULL;

with, and change the lines

CombineMatchingInput combined;
matching = combined.getHook(pythia);

to

matching = new JetMatchingMadgraph();

This will allow you to access the Pythia-internal DJRs through

vector<double> my_DJRs = matching->getDJR();

after the pythia.next() call (note that the function is called "getDJR" in later Pythia8.2 versions, and "GetDJR" in earlier versions, so that you might have to adjust, or update your Pythia installation). Furthermore, you can retrieve the "number of additional hard partons" via

int my_nHardPartons = matching->nMEpartons().first;

This should allow you to histogram the DJRs, split into jet bins. In case you need inspiration on how to implement the histograms, please have a look at

http://slac.stanford.edu/~prestel/Support/main89-djrtest.cc

Hope that helps,
Stefan

Jack Y. Araz (jackaraz) said : #9

Dear Stefan, Oliver

Thank you very much, probably the problem is with ma5 because DJR output from py8 looks as expected. I have further question about x-section difference between py8 and mg5 if you don't mind.

Since I'm getting drastic difference between mg5-py6-py8, wanted to ask if its normal and which one should I rely on. In order to get fast response I generated same event as I mentioned above with just neutrinos and the result that I get for 100 events is

mg5 x-sec: 2.9e+04 pb
py6 x-sec: 9507 pb
py8 x-sec: 6347 pb (this result was around 7484 pb for 10000 events but other values were similar)

so to understand this difference I checked PDF set which is nn23lo1 for mg5 thus I set PDF:pSet = 13. Also I fixed renormalization and factorization scale in mg5 ( T = fixed_ren_scale , T = fixed_fac_scale) which is also done in py8 (TimeShower:alphaSvalue = 0.11798 (I also tried 0.130 but changed since mg5 uses alphaS value provided by PDF set as far as I know), TimeShower:alphaSorder = 0, TimeShower:alphaEMorder = 0, same for space shower) also maximum jet flavor set to 5 (5 = maxjetflavor). The only thing that comes to my mind is quark masses which generally taken as zero but I'm taking them nonzero both in py8 and mg5 thus since they are identical it doesn't make sense to me to get such difference

1:m0 = 0.0048
2:m0 = 0.0023
3:m0 = 0.095
4:m0 = 1.275
5:m0 = 4.66
6:m0 = 173.21
11:m0 = 0.000510999
13:m0 = 0.105658
15:m0 = 1.77686

Thus I couldn't be sure if i'm making a mistake along the line or which x-sec value should I rely on.

Thank you very much
Best regards
Jack

Hi,

Not sure what you are comparing.
The cross-section in mg5 is likely to be unphysical since it still has the double counting between the various multiplicities.
The correct cross-section is only available after the shower when the correct jet veto has been applied to remove this double counting.

Cheers,

Olivier

Jack Y. Araz (jackaraz) said : #11

Hi Oliver,

I'm just trying to understand the x-sec difference in py6-py8 after shower and correct jet veto, or in other words I'm trying to set my self a benchmark to understand if I'm doing something wrong while generating the bkg. Where the results are quite different even after checking the key points that you mentioned in other questions. Thus I was expecting to get somewhat closer result from py6 and py8, I believe %33 difference is quite a lot.

Thank you very much
Best regards
Jack

Stefan Prestel (prestel) said : #12

Dear Jack,

Before being able to say anything, I would need more information. Are you positive that your xqcut and JetMatching:qCut have identical values in your two simulations? Are you using Pyitha6 with Q^2- or pT-ordered showers? Did you check that Pythia8 does not try to read more events than available from the MG5-produced LHE file? To answer these questions, it would be helpful if you could provide your run_card.dat, the first ~1000 lines of the LHEF, your Pythia 8 input card, and the Pythia 8 terminal output. That should help us figure out if something is amiss.

Cheers,
Stefan

Jack Y. Araz (jackaraz) said : #13

Dear Stefan,

> Are you positive that your xqcut and JetMatching:qCut have identical values in your two simulations?

Yes, they are.

>Are you using Pyitha6 with Q^2- or pT-ordered showers?

I'm using default pythia_card as done in TASI lectures 2013

>Did you check that Pythia8 does not try to read more events than available from the MG5-produced LHE file?

Yes, I'm initiating the number.

>To answer these questions, it would be helpful if you could provide your run_card.dat, the first ~1000 lines of the LHEF, your Pythia 8 input card, and the Pythia 8 terminal output. That should help us figure out if something is amiss.

You can find the files that I used in this link: https://www.dropbox.com/sh/0xujglzxp1t7cki/AADvOHwyUGVxKvrEMS0wbek6a?dl=0

I'm constantly changing the pythia8 card but you can find the general frame that I used and one of the terminal output is there. I also included the code that I used with your advice. Hope it helps, I really want to understand where that I'm doing it wrong.

Thank you very much
Best regards
Jack

Stefan Prestel (prestel) said : #14

Hi Jack,

From your inputs, it seems that

- you're requesting 60K events from Pythia, whereas only 20K are stored in the LHEF. The normalization will thus be off. Maybe try requesting 10K events.
- Since I'm not sure if (or how) Pythia 6 sets etaJetMax, please set etaJetMax = 100
- I'm fairly certain that Pythia 6 will use qcut = 1.5*xqcut, although I can't quickly find the TASI 2013 setup on the web. So please try setting qCut = 15, which is (for MLM jet matching) a sensible choice if xqcut = 10.

Then some further comments:

- Your value of exclusive matches the default. You thus do not need to put this setting explicitly.
- With the exclusive value you use, you don't need to set nJet

I hope this helps.

Best
Stefan

Jack Y. Araz (jackaraz) said : #15

Dear Stefan

Thank you very much, I didn't realised that nEvents has an effect on x-section. it also changed the final result when I combine all boson jet events.

Thanks for your time
Best regards
Jack

Jack Y. Araz (jackaraz) said : #16

Dear Stefan, Oliver

I want to ask the same question for a supersymmetric model this time. I want to do the DJR analysis for a U(1)' extended MSSM which created by SARAH-SPHENO (note that spheno doesnt support block qnumbers can the problem related to that?). Basically in a similar manner that I mentioned above I generated

generate p p > DM DM
add process p p> DM DM j
add process p p > DM DM j j

for xqcut = 10, I generated 10^8 event in MG5 with a x-sec : 0.19 pb (increases with #events) Then I gave it to py8 as done above (files can be found from this link: https://www.dropbox.com/sh/0xujglzxp1t7cki/AADvOHwyUGVxKvrEMS0wbek6a?dl=0) However I can not see any DJR plots due to very low x-sec which is around 10^-40. So I also tried without jet matching but this time I couldn't generate more than 200 events. I tried with higher and lower nEvent values but low values doesnt generate enough event and high values start to abort due to somany error such as: 1393 Error in Pythia::check: unknown particle code

Jack Y. Araz (jackaraz) said : #17

sorry the rest didnt appeared:
...
Is there any function that I can introduce a particle or should I do in in entirely different way?
Thank you very much
Best regards
Jack

Can you help with this problem?

Provide an answer of your own, or ask Jack Y. Araz for more information if necessary.

To post a message you must log in.