[Rivet] Pythia+PU+Rivet problem

[Andy Buckley]

Hi again,

Thanks for the extra info. Since the exception is being thrown in the
AnalysisHandler::analyze, but before the Analysis::analyze, I think it's
very likely that it's coming from this checking code:

if (!compatible(beams, _beams) || !fuzzyEquals(sqrts, sqrtS())) {
  cerr << "Event beams mismatch: " << PID::toBeamsString(beams)
       << " @ " << sqrts/GeV << " GeV" << " vs. first beams "
       << this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl;
  exit(1);
}

specifically the PID::toBeamsString(beams) call. My suspicion is that
this is because PYTHIA6 doesn't properly identify the beam particles,
and the convention that (if not explicitly declared) they are the first
two particles in the event is being broken by the pile-up simulation.

This seems particularly likely because I'm fairly sure that event
printout is from PYTHIA6, not from Rivet or HepMC (the "K" headers look
familiar). So I think PYTHIA itself is recognising that the pile-up mode
is producing problem events, but not vetoing them...

Andy


On 15/04/14 10:42, Tomasz Fruboes wrote:
> Hi,
> 
> 
>  please find part of the printout attached below. After inserting
> additionall printouts (and protecting all of my printouts by
> std::cout.flush()) I see, that the  exception occurs after the
> Rivet::AnalysisHandler::analyze() call and before calling the analyze()
> method of my analysis class.
> 
>  I have checked - the sumOfWeight counter is not to affected by this (ie
> it's lower than the number of generated events by number of exceptions).
> 
>  Cheers,
>    Tomasz
> 
> 
> 
> 
> 
>                             Event listing (no momenta)
> 
>     I  particle/jet  K(I,1)   K(I,2) K(I,3)     K(I,4)      K(I,5)    C
> tag  AC tag
> 
>     1  !p+!              21      2212    0           0           0   
> 0       0
>     2  !p+!              21      2212    0           0           0   
> 0       0
>     3  !u!               21         2    1           0           0   
> 1       0
>     4  !ubar!            21        -2    2           0           0   
> 0       3
>     5  !u!               21         2    3           0           0   
> 0       0
>     6  !ubar!            21        -2    4           0           0   
> 0       0
>     7  !g!               21        21    0           0           0   
> 0       0
>     8  !Z0!              21        23    0           0           0   
> 0       0
>     9  !mu-!             21        13    8           0           0   
> 0       0
>    10  !mu+!             21       -13    8           0           0   
> 0       0
>    11  (Z0)              11        23    8          12          13   
> 0       0
>    12  mu-                1        13    9           0           0   
> 0       0
>    13  mu+                1       -13   10           0           0   
> 0       0
>    14  (ud_1)        A   12      2103    1          63          63   
> 0       0
>    15  (g)           I   12        21    1          63          63   
> 0       0
>    16  (g)           I   12        21    2          63          63   
> 0       0
> 
> 
> On 15.04.2014 10:20, Andy Buckley wrote:
>> On 14/04/14 23:17, Tomasz Fruboes wrote:
>>> Hi,
>>>
>>>   thanks for answer. No, the analysis code doesnt access the particle
>>> names. I'm getting quite a lot event printouts (not requested :) ), so
>>> I'm guessing that rivet is finding something weird with the events and
>>> prints it out as a debug info. And time to time, during the printout, he
>>> asks for the name of pdg=1...3 resulting with the crash.
>>
>> Sounds reasonable, but I don't know where in Rivet there would be an
>> event printout... can you send a bit of the printout so we can tell if
>> it's coming from Rivet or HepMC (or other) code? Strange...
>>
>> If you put printouts in your code (or some other method), can you tell
>> if the exception is being thrown before your code executes, or if it's
>> triggered by something that you are calling in the analysis?
>>
>>>   PU in Pythia6 seems to be easiest to get on generator level (AFAIK in
>>> pythia8 this has to be handled on c++ level - by defining two separate
>>> pythia sources - and cannot be controlled from "cards" level). BTW the
>>> workaround is not intended to be included in any official release (I'm
>>> not a developer/maintainer of the rivet wrapper present in CMSSW), this
>>> is just to keep going with work. So - do you think this can be dangerous
>>> on the rivet side (ie can broke some rivet internals)?
>>>
>>>   My main concern with the workaround is that event weight counting can
>>> get corrupted. I'll try to investigate this (unless you got immediate
>>> answer :) )
>>
>> I think that is possible, yes. And depending on where the error occurs,
>> you might be filling some histograms but not others.
>>
>> Andy
>>
>>
>>
>>> On 14.04.2014 22:13, Andy Buckley wrote:
>>>> Hi Tomasz,
>>>>
>>>> I didn't know anyone was using the old PYTHIA pile-up system!
>>>> Looking in
>>>> the Rivet code, that exception is thrown when it tries to find a text
>>>> name for an unknown particle: since PIDs 1,2,3 are all quarks and hence
>>>> should never be seen in a stable event final state, it leads to your
>>>> error.
>>>>
>>>> Do you know why your code could be trying to access Particles' string
>>>> names? It is possible that there's a bug in PYTHIA6 when run this way,
>>>> which leads to partons appearing in the final state: I'm not sure.
>>>> Obviously it's better not to put catch-all exception handling around
>>>> the
>>>> analyze function call -- you can hide serious problems that way and
>>>> make
>>>> debugging hard for everyone: presumably the error is being thrown by
>>>> something inside your analysis routine, so maybe wrap *that* in the try
>>>> block (and make the try block as small as possible) rather than the
>>>> general CMSSW Rivet wrapper?
>>>>
>>>> By the way, I strongly recommend updating to use Rivet 1.9.0 -- your
>>>> analysis should continue to work as before, but several bugs were fixed
>>>> between 1.8.2 and that. All future releases, bugfixes, etc. are in the
>>>> Rivet 2.x series.
>>>>
>>>> Best wishes,
>>>> Andy
>>>>
>>>>
>>>> On 14/04/14 16:52, Tomasz Fruboes wrote:
>>>>> Dear Experts,
>>>>>
>>>>>    in our rivet (1.82) analysys I'm trying to include pileup.
>>>>> Pileup is
>>>>> added on generator level by pythia6 using following settings:
>>>>>
>>>>>               'MSTP(131)=1
>>>>>               'MSTP(132)=3
>>>>>               'MSTP(133)=1
>>>>>               'PARP(131)=0.15
>>>>>
>>>>> Above gives n=10 PU interactions on average. In every 1 out of ~300
>>>>> events I'm getting following exception
>>>>>
>>>>> A std::exception was thrown.
>>>>> Particle ID '1' not known.
>>>>>
>>>>> (where ID can be 2 or 3 as well).
>>>>>
>>>>> I have added try/catch block in order to handle this (see [1] below).
>>>>> Could you comment if this is a safe workaround? Or maybe there is a
>>>>> better solution?
>>>>>
>>>>>    Thank you,
>>>>>     Tomasz Fruboes
>>>>>
>>>>>
>>>>> [1] Source code of our rivet wrapper:
>>>>>
>>>>> https://github.com/cms-sw/cmssw/blob/CMSSW_7_1_X/GeneratorInterface/RivetInterface/plugins/RivetAnalyzer.cc
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> Line ~120 was changed to
>>>>>     //run the analysis
>>>>>     try {
>>>>>         _analysisHandler.analyze(*myGenEvent);
>>>>>     } catch (...) {
>>>>>           std::cout << "Exception caught...\n" ;
>>>>>     }
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> Rivet mailing list
>>>>> Rivet at projects.hepforge.org
>>>>> https://www.hepforge.org/lists/listinfo/rivet
>>>>
>>>>
>>>
>>
>>
> 


-- 
Dr Andy Buckley, Royal Society University Research Fellow
Particle Physics Expt Group, University of Glasgow / PH Dept, CERN