Jet matching/merging & DJR analysis with pythia8
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:/
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:
JetMatching:
JetMatching:nJetMax = 4
JetMatching:
Check:epTolErr = 1e-2
LHEFInputs:nSubruns = 2
Main:subrun = 0
Beams:LHEF = <l+l- events PATH>/unweighte
Main:subrun = 1
Beams:newLHEFsa
Main:LHEFskipInit = on
Beams:LHEF = <vl vl~ event PATH>/unweighte
in ma5:
./bin/ma5 -H
submset main.merging.check = true
submset main.merging.njets = 4
import <PATH>/
submit
Note that I also tried to use 'set main.merging.
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
- Assignee:
- Stefan Prestel Edit question
- Last query:
- 2016-08-11
- Last reply:
- 2016-07-30
This question was reopened
- 2016-08-11 by Jack Y. Araz
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:/
>
> 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,
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:/
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
...
JetMatchingMadgraph *fJetMatchingPy
fJetMatchingPy8
pythia.
const std::vector<double> djrmatch = fJetMatchingPy8
std::ofstream myoutput;
myoutput.
for(unsigned int i=0; i < std::min(
{
myoutput << djrmatch[i]<< "\t"
<< fJetMatchingPy8
<< fJetMatchingPy8
<< std::endl;
}
where JetMatchingPy8I
main89_
const std::vector<double> djrmatch = fJetMatchingPy8
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
JetMatchingMadg
with, and change the lines
CombineMatching
matching = combined.
to
matching = new JetMatchingMadg
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-
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://
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:
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:/
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:/
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.