Don Benton pointed out at Ann Arbor that if you run the jet finder on
the MC truth as output by Gismo, using the code given as an example at
http://www-sldnt.slac.stanford.edu/jas/documentation/lcd/tutorial9/defau
lt.htm
it produces some jets which have unphysically high energies. This is
caused by a problem in the Gismo MC truth table, where it lists the same
electron many times, and in each case classifies it as INTERACTED
(==FINALSTATE). The problem can be seen very clearly by running the
JetFinder on Bhabha events.
A real solution needs a fix in gismo, but as a workaround it is possible
to test not only that the particle is listed as FINALSTATE, but also
that it has no daughters. Thus we take only the final electron in the
list as the "finalstate" electron. A simple job to show how to do this
is attached.
Note that the FastMC built into hep.lcd uses the same test as the
original JetFinder for "FinalState" particle, and will thus exhibit a
similar problem if run on Gismo files (as opposed to StdHEP files).
Tony
import hep.analysis.*;
import hep.lcd.event.*;
import hep.lcd.util.predicate.*;
import hep.lcd.util.driver.*;
import hep.physics.*;
import hep.lcd.physics.*;
final public class JetDriver extends Driver
{
public JetDriver()
{
add(new JetFinder());
out.println("Remote server running hep.lcd version "+
hep.lcd.Version.getVersion());
}
}
final class JetFinder extends AbstractProcessor
{
IJetFinder jf = new DurhamJetFinder(0.02);
Predicate predicate = new Predicate()
{
public boolean accept(Object o)
{
Particle p = (Particle) o;
if (p.getStatusCode() != Particle.FINALSTATE)
return false;
if (p.getDaughters().hasMoreParticles()) return
false;
return true;
}
};
public void process(final LCDEvent event)
{
// Fill histograms and/or scatter plots here.
ParticleVector pv = event.getMCParticles();
Filter f = new Filter(pv.particles(),predicate);
jf.setEvent(f);
histogram("njets").fill(jf.njets());
histogram("nparts").fill(pv.size());
int np = 0;
double etot = 0;
for (int i=0; i<jf.njets(); i++)
{
np += jf.nParticlesPerJet(i);
etot += jf.jet(i).t();
}
histogram("np").fill(np);
histogram("etot").fill(etot);
}
}
|