Jet matching/merging & DJR analysis with pythia8

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

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.

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

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

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

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

-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

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

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

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

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

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

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

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

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