package hep.physics.jet;
import hep.physics.filter.Predicate;
import hep.physics.vec.BasicHepLorentzVector;
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.HepLorentzVector;
import hep.physics.vec.VecOp;
import java.util.List;
import java.util.ArrayList;
import java.util.Collection;
import java.util.Iterator;
/**
* Base class for jet finders
* @author G.Bower
* @version 12/6/98
*/
public abstract class AbstractJetFinder implements JetFinder
{
private double defaultMassSquared = 0;
private final static int UNASSOC = -999;
private int m_injets;
// m_jet[i] is the 4 vector sum of all the particles in jet i.
private HepLorentzVector[] m_jet;
//m_ipart_jet_assoc[i] = j means particle i was placed in jet j.
private int[] m_ipart_jet_assoc;
// m_inparts_per_jet[i] = j means jet i has j particles in it.
private int[] m_inparts_per_jet;
// m_ifewest_tracks is number of tracks in the jet with the fewest tracks.
private int m_ifewest_tracks;
//
private double m_dycut;
private boolean m_resultsValid = false;
private List m_in;
private List m_4vec;
private double m_devis = 0;
//
protected int m_np;
protected HepLorentzVector[] m_part;
//
protected double[][] ymass;
//
abstract double masscut(double ycut, double evis, double esum);
//
abstract double calculate_mass(HepLorentzVector p1, HepLorentzVector p2);
//
abstract void combine_particles(int im, int jm);
//
public int njets()
{
if (!m_resultsValid) doFindJets();
return m_injets;
}
public HepLorentzVector jet(int index)
{
if (!m_resultsValid) doFindJets();
return m_jet[index];
}
public List particlesInJet(int index)
{
if (!m_resultsValid) doFindJets();
List result = new ArrayList();
for (int i=0; i masscut(m_dycut,m_devis,esum)) break;
//
// combine particles im and jm update associations
//
combine_particles(im,jm);
m_ipart_jet_assoc[jm] = im;
for(int j = 0 ; j < m_np ; j++ )
{
if( m_ipart_jet_assoc[j] == jm )
{
m_ipart_jet_assoc[j] = im;
}
}
//
// Recalculate a mass for all pairs that contain im and a remaining
// UNASSOC particle(jet). Also recalculate esum.
//
esum = 0.0;
for(int j = 0 ; j < m_np ; j++ )
{
if (m_ipart_jet_assoc[j] != UNASSOC ) continue;
esum = esum + m_part[j].t();
if( j == im) continue;
int imin = Math.min(j,im);
int imax = Math.max(j,im);
double cmass = calculate_mass(m_part[imin], m_part[imax]);
if ( cmass >= 0.0 )
{
ymass[imin][imax] = cmass;
}
}
}
//
// finish up by filling jet array.
//
for(int i = 0 ; i < m_np ; i++ )
{
if (m_ipart_jet_assoc[i] == UNASSOC) m_injets++;;
}
m_jet = new HepLorentzVector[m_injets];
m_inparts_per_jet = new int[m_injets];
int nj = 0;
int ntrk;
m_ifewest_tracks = Integer.MAX_VALUE;
for(int i = 0 ; i < m_np ; i++ )
{
if (m_ipart_jet_assoc[i] != UNASSOC) continue;
m_jet[nj] = m_part[i];
ntrk = 1;
for (int j = 0 ; j < m_np ; j++)
{
if(m_ipart_jet_assoc[j] == i)
{
m_ipart_jet_assoc[j] = nj;
ntrk++;
}
}
m_ipart_jet_assoc[i] = nj;
m_inparts_per_jet[nj] = ntrk;
if( ntrk < m_ifewest_tracks) m_ifewest_tracks = ntrk;
nj++;
}
}
/**There are a variety of collinear, infra-red safe jet finders available.
* The differences between them are
* (1) How they determine the cutoff mass from the ycut parameter.
* (2) How they combine two four-vectors together as they build up the jet.
* (3) what mass they assign to a four-vector.
* The following methods are the algorithms used in these jet finders.
* A concrete jet finder is an extension of the AbstractJetFinder (the current
* class) that specifies a masscut method, a combine_particles method and
* a calculate_mass method.
*/
//
// Cutoff mass (masscut) method.
//
protected double standard_masscut(double ycut, double evis)
{
return ycut*evis*evis;
}
protected double geneva_masscut(double ycut)
{
return ycut;
}
protected double jadeP0_masscut(double ycut, double esum)
{
return ycut*esum*esum;
}
//
// combine_particles methods
//
protected void four_vector_combine( int im, int jm)
{
m_part[im] = VecOp.add( m_part[im], m_part[jm] );
}
protected void jadeP_combine( int im, int jm )
{
Hep3Vector v = VecOp.add( m_part[im].v3(), m_part[jm].v3() );
m_part[im] = new BasicHepLorentzVector( v.magnitude(), v);
}
protected void jadeE0_combine( int im, int jm )
{
HepLorentzVector v = VecOp.add( m_part[im], m_part[jm] );
double ekeratio = v.t()/v.v3().magnitude();
m_part[im] = new BasicHepLorentzVector(v.t(),VecOp.mult(ekeratio, m_part[im].v3()));
}
//
// calculate_mass methods
//
protected double four_vector_mass(HepLorentzVector part1, HepLorentzVector part2)
{
return VecOp.add(part1, part2).magnitudeSquared();
}
protected double jade_mass(HepLorentzVector part1, HepLorentzVector part2)
{
double e1 = part1.t();
double e2 = part2.t();
Hep3Vector v1 = part1.v3();
Hep3Vector v2 = part2.v3();
double costh = VecOp.dot(v1,v2)/(v1.magnitude()*v2.magnitude());
return 2*e1*e2*( 1 - costh );
}
protected double durham_mass(HepLorentzVector part1, HepLorentzVector part2)
{
double e1 = part1.t();
double e2 = part2.t();
Hep3Vector v1 = part1.v3();
Hep3Vector v2 = part2.v3();
double costh = VecOp.dot(v1,v2)/(v1.magnitude()*v2.magnitude());
double lessorE = Math.min(e1,e2);
return 2 * lessorE*lessorE * ( 1 - costh );
}
protected double geneva_mass(HepLorentzVector part1, HepLorentzVector part2)
{
double e1 = part1.t();
double e2 = part2.t();
if ( e1 == 0.0 & e2 == 0.0 ) return -9999.0;
Hep3Vector v1 = part1.v3();
Hep3Vector v2 = part2.v3();
double costh = VecOp.dot(v1,v2)/(v1.magnitude()*v2.magnitude());
return (8.0/9.0)*(1 - costh)*(e1*e2)/((e1 + e2)*(e1 + e2));
}
}
Ads help maintain this website.