HepMC from Herwig++

[holdom]

I use Herwig++ 2.6.3 to generate HepMC events with W+W- yielding e and mu.
When parton showers, hadronization and decay are all turned on then DelphesHepMC gives
ERROR: invalid vertex format
When decay is turned off then DelphesHepMC succeeds.
When hadronization is turned off as well then muons no longer appear in DelphesHepMC output.
When parton showers are turned off as well then the muons reappear.

[Pavel Demin]

In the previous ticket #125, it was concluded that it's Herwig++ that output numbers outside of the double precision floating point number range (~2.2e−308, ~1.8e308).

In the file that was attached to that ticket I found the following line

V -237 0 2.8878950987394353e-314 1.1789752640139257e-313 2.2493714921184413e-313 0 0 1 0

e-314 is clearly outside the range.

I would suggest reporting this problem to Herwig++ developers.

For Delphes 3, I've found a solution that returns the same error on MacOSX and on Linux when it finds the illegal double precision floating point numbers:

** ERROR: invalid vertex format

** Exiting...

In the same ticket #125, a workaround was suggested (to set to 0.0 everything that is less than 1.0e-300):

1. Put the following lines in a perl script 'script.pl':

   @newline = ();
   foreach $token (@F) {
     if($token =~ /e-/ && abs($token) < 1e-300) {
       push(@newline, "0.0e+0");
     }
     else {
       push(@newline, $token);
     }
   }
   print "@newline\n";
   

2. Then invoke this script with the following command:

   perl -lan script.pl hepmc_input.hepmc | ./DelphesHepMC examples/delphes_card_CMS.tcl delphes_output.root
   

[Pavel Demin]

Concerning the appearing and disappearing muons, I suspect that it's the status of the particles that changes. By default, I suppose that final state particles have status 1.

I could find the following information in the Herwig++ documentation:

​https://herwig.hepforge.org/trac/wiki/FaQs#WhichHepMCstatuscodesareusedWhatdotheymean

When hadronization is turned off, muons are probably considered as intermediate particles and get status code 2 or 11-200.

Could you, please check the muons status codes in your HepMC files?

[Pavel Demin]

Another workaround for the out-of-range numbers could be replacing

  return start != fBuffer && errno != ERANGE;

with

  return start != fBuffer;

on the lines 38 and 48 in classes/DelphesStream.cc

After this modification, the too small values will be replaced with 0 and too large values will be replaced with HUGE_VAL or -HUGE_VAL.

[Pavel Demin]

In any case, this problem should be reported to the Herwig++ developers at

​https://herwig.hepforge.org/trac/

It's not normal that Herwig++ outputs unphysical numbers.

[Pavel Demin]

Just found a warning about unphysical values in the Herwig++ HepMC output:

​https://herwig.hepforge.org/trac/wiki/FaQs#HowdoImakeHerwigwriteouteventsasHepMC

So, they won't probably fix this problem.

I've modified classes/DelphesStream.h and classes/DelphesStream.cc to output a warning the first time it encounters an out-of-range value.

Could you, please check whether new version works for you

​http://cp3.irmp.ucl.ac.be/downloads/Delphes-3.0.4.tar.gz

[Pavel Demin]

I've just realized that the ticket notification system did not work.

The notifications are fixed now.

However, you were not notified about new comments to this ticket.

Could you, please check the comments 3 and 5 via the Delphes web site?

[holdom]

Yes, Delphes-3.0.4 does now succeed on hadronized and decayed events from Herwig++, and I see the warning.

[Pavel Demin]

Thank you for the test!

[holdom]

As for the muons, the non-hadronized HepMC file seems to have muons with status codes 1 and 11, just as for the electrons (assuming I know how to read the file). But muons don't show up in the LHCO file (or the root file as far as I can tell--I usually only work with LHCO files). The same HepMC file fed through Delphes 2 does produce muons in the LHCO file, so something has changed.

[Pavel Demin]

Could you, please send me a link to one of your non-hadronized HepMC files?

I'll try to understand what is going on.

[holdom]

I attached the file.

[Pavel Demin]

Thank you for the file.

The production vertices of all muons with status 1 are far outside of the detector volume (defined as cylinder with radius ~1m and diameter ~3m).

Here are the first 2 muons with status 1 from the file (with the V line describing the production vertex position):

V -28 0 3.2906898720063623e+07 1.3936505537842816e+08 -9.4998864093059506e+06 1.4351225602307361e+08 0 1 0
P 10042 -13 1.8040581294913210e+01 7.6404240722072984e+01 -5.2081320247608600e+00 7.8677864590809861e+01 1.0565836900000000e-01 1 0 0 0 0

V -26 0 1.8446232446488759e+08 -4.1742524862695321e+07 1.0790521272128012e+09 1.0955009906088841e+09 0 1 0
P 10040 -13 5.2153744216572264e+01 -1.1802024998645134e+01 3.0508456836516001e+02 3.0973521893405064e+02 1.0565836900000000e-01 1 0 0 0 0

Then we have the following expression in the ParticlePropagator module (modules/ParticlePropagator.cc):

    // check that particle position is inside the cylinder
    if(TMath::Hypot(x, y) > fRadius || TMath::Abs(z) > fHalfLength)
    {
      continue;
    }

It means that the particles originating from a vertex outside of the detector volume don't participate in the simulation that follows the particle propagator.

[holdom]

Okay that may be the case for those particular muons, but Delphes 2 sees the muons from W decay while Delphes 3 does not. I would rather have the Delphes 2 behavior.

[Pavel Demin]

The problem is that in the file that you provide all the muons with status 1 have similar characteristics and have their production vertices far (kilometers) outside of the detector.

I've just checked how particle propagation is done in Delphes 2 (src/BFieldProp.cc).

Looks like Delphes 2 propagates all the particles that are produced outside of the detector volume with a straight line drawn from (0, 0, 0) without taking into account the coordinates of the production vertex.

Here is the code from src/BFieldProp.cc:

void TrackPropagation::bfield(TRootGenParticle *Part) {

  // initialisation, valid for z_max==0, R_max==0 and q==0
  Part->EtaCalo = Part->Eta;
  Part->PhiCalo = Part->Phi;//-atan2(Part->Px,Part->Py);

// skipped a few lines

  // out of tracking coverage?
  if(sqrt(X*X+Y*Y) > R_max){return;}
  if(fabs(Z) > z_max){return;}

It sets the propagated values of eta and phi (EtaCalo and PhiCalo) before checking if the particle is within the tracker volume.

Then EtaCalo is used to check if the particle is within the tracker volume => looks like a bug.

Here is the corresponding code from Delphes.cpp:

              if( particle->Charge !=0 &&
                  fabs(particle->EtaCalo)< DET->CEN_max_tracker && // stays in the tracker -> track available
                  ( DET->FLAG_bfield ||
                    (!DET->FLAG_bfield && particle->PT > DET->TRACK_ptmin)
                    )
                  )
                {
                  bool trackRec=true;
                  if((grandom->Uniform()*100.) > DET->TRACK_eff)trackRec=false;
                  // 2.1a.2.3.1 Filling the particle properties + smearing
                  // Hypothesis: the final eta/phi are the ones from the generator, thanks to the track reconstruction
                  // This is the EnergyFlow hypothesis
                  particle->SetEtaPhi(particle->Eta,particle->Phi);
                  float sET=UNDEFINED; // smeared ET, computed from the smeared E -> needed for the tracks

                  // 2.1a.2.3.2 Muons
                  if (pid == pMU && fabs(particle->EtaCalo)< DET->CEN_max_mu && trackRec==true)
                    {
...

I'm not sure if we should call this behavior correct.

[holdom]

I don't understand because if I turn on hadronization in Herwig++ then Delphes 3 sees the muons coming from the decays of W+W- that are produced in each event. In each event one W decays to mu and the other decays to e. Why should these muons disappear when I turn off hadronization?

[Pavel Demin]

Could you, please check the coordinates of the production vertices for the muons with status 1 when hadronization is on?

[holdom]

I have attached the hadronized file. I have trouble looking at such files.

[Pavel Demin]

Thank you for the hadronized file.

The coordinates of the production vertices for muons are within the detector volume:

V -33 0 5.6718499319230702e-15 1.0131239098778490e-12 -5.5129365158446199e-13 1.4166253965734495e-12 0 1 0
P 10052 -13 1.8040581294913210e+01 7.6404240722072984e+01 -5.2081320247608600e+00 7.8677864590809861e+01 1.0565836900000000e-01 1 0 0 0 0

In the HepMC file without hadronization the muon production loooks as

W+ -> gamma nu mu+ with status code 11 -> gamma nu mu+ with status code 1

For gamma and nu, vertex position is OK. For mu it's ~1e7 mm.

Is it possible that Herwig++ for some reason considers mu as unstable and tries to decay it?

Herwig++ documentation mentions some parameters that define stable particles
​https://herwig.hepforge.org/trac/wiki/FaQs#HowdoImakeagroupofparticlesstable

I've seen that some Tevatron and LHC configuration files for Herwig++ contain

set /Herwig/Decays/DecayHandler:LifeTimeOption 0
set /Herwig/Decays/DecayHandler:MaxLifeTime 10*mm

Do you have something similar in you Herwig++ configuration?

[holdom]

I have attached a file where hadronization is turned off, the decay handler is turned on, and I include your two decay handler statements above. Delphes 3 still sees no muons. They again reappear when hadronization is turned on.

[Pavel Demin]

In the PHd.hepmc file all muons with status 1 have their production outside the detector volume:

V -28 0 3.2906898720063623e+07 1.3936505537842816e+08 -9.4998864093059506e+06 1.4351225602307361e+08 0 1 0
P 10042 -13 1.8040581294913210e+01 7.6404240722072984e+01 -5.2081320247608600e+00 7.8677864590809861e+01 1.0565836900000000e-01 1 0 0 0 0

[Pavel Demin]

I'm out of ideas.

Here is the current status of this ticket:

1. Delphes 3 is modified to accept out of range values in the HepMC files. Too small values are reset to 0 and too large positive and negative values are reset to inf and -inf.

2. Delphes 2 erroneously propagates all the particles that are produced outside of the detector volume with a straight line drawn from (0, 0, 0). Delphes 3 rejects such particles.

3. When hadronization is turned on, Herwig++ output files can be processed with Delphes 3.

4. When hadronization is turned off, Herwig++ output files contain muons produced at ~1e4 meters from the interaction point. Delphes 3 reject these muons.

I would propose to address the following questions to the Herwig++ developers:

1. How to configure Herwig++ to output HepMC files with valid floating point numbers?

2. How to configure Herwig++ to output muons produced at the interaction point when hadronization is turned off?

[holdom]

Yes, it would be great if Herwig++ developers knew about these problems.

In the meantime if I was to modify my local Delphes code in modules/ParticlePropagator.cc that you mentioned above so as to accept particles outside the detector volume, I would presumably recover the muons. Would this work? Are there unwanted side effects?

[Pavel Demin]

Yes, it would be great if Herwig++ developers knew about these problems.

Could you, please contact them?

In the meantime if I was to modify my local Delphes code in modules/ParticlePropagator.cc that you mentioned above so as to accept particles outside the detector volume, I would presumably recover the muons.
Would this work? Are there unwanted side effects?

To recover the muons and to reduce eventual unwanted side effects, you can replace lines 122-125 in modules/ParticlePropagator.cc:

    if(TMath::Hypot(x, y) > fRadius || TMath::Abs(z) > fHalfLength)
    {
      continue;
    }

with

    if(TMath::Hypot(x, y) > fRadius || TMath::Abs(z) > fHalfLength)
    {
      if(TMath::Abs(candidate->PID) == 13)
      {
        x = 0.0;
        y = 0.0;       
        z = 0.0;               
      }
      else
      {
        continue;
      }
    }