lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -N HelixUtils.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HelixUtils.java 19 Jun 2009 19:11:19 -0000 1.1
@@ -0,0 +1,430 @@
+/*
+ * HelixUtils.java
+ *
+ * Created on February 1, 2008, 5:48 PM
+ *
+ */
+
+package org.lcsim.contrib.Mbussonn;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.matrix.BasicMatrix;
+import hep.physics.matrix.Matrix;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import org.lcsim.fit.circle.CircleFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.fit.helicaltrack.TrackDirection;
+
+/**
+ * Assorted helix utilities that operate on HelicalTrackFits.
+ *
+ * Since these utilities are typically called many times for a
+ * given HelicalTrackFit (or in some cases for a given CircleFit),
+ * a number of derived quantities used by the methods in this class
+ * are cached to minimize re-computation. Since a HelicalTrackFit
+ * (or CircleFit) is immutable, these quantities only need to be
+ * calculated once.
+ * @author Richard Partridge
+ * @version 1.0
+ */
+public class HelixUtils {
+ private static double _minslope = 1.e-6;
+ private static double _maxpath = 2400;
+ private static double _epsilon = 0.01; // value below witch a z vector Or an angle is considered as null
+
+ /**
+ * Creates a new instance of HelixUtils.
+ */
+ public HelixUtils() {
+ }
+
+ /**
+ * Return the x-y path length between two HelicalTrackHits.
+ * @param cfit CircleFit to be used in calculating the path length
+ * @param hit1 first hit
+ * @param hit2 second hit
+ * @return path length from hit1 to hit2
+ */
+ public static double PathLength(CircleFit cfit, HelicalTrackHit hit1, HelicalTrackHit hit2) {
+ // Return the path length between hit1 and hit2
+ return PathCalc(xc(cfit), yc(cfit), RC(cfit), hit1.x(), hit1.y(), hit2.x(), hit2.y());
+ }
+
+ /**
+ * Return the x-y path length from the DCA.
+ * @param cfit CircleFit to be used in calculating the path length
+ * @param hit hit to be used for the path length calculation
+ * @return path length from the DCA to the hit
+ */
+ public static double PathLength(CircleFit cfit, HelicalTrackHit hit) {
+ // Return the path length from the DCA
+ return PathCalc(xc(cfit), yc(cfit), RC(cfit), x0(cfit), y0(cfit), hit.x(), hit.y());
+ }
+
+ /**
+ * Return the x-y path length from the DCA to a HelicalTrackHit.
+ * @param helix HelicalTrackFit to be used in calculating the path length
+ * @param hit hit to be used for the path length calculation
+ * @return path length from the DCA to the hit
+ */
+ public static double PathLength(HelicalTrackFit helix, HelicalTrackHit hit) {
+ // Return the path length from the DCA
+ return PathCalc(helix.xc(), helix.yc(), helix.R(), helix.x0(), helix.y0(), hit.x(), hit.y());
+ }
+
+ /**
+ * Return the x-y path length to a z-plane.
+ * @param helix HelicalTrackFit to be used in calculating the path length
+ * @param z location of z-plane
+ * @return path length from the DCA to the z-plane
+ */
+ public static double PathToZPlane(HelicalTrackFit helix, double z) {
+ // Find the z distace between the DCA and the z plane and calculate the x-y path length
+ double zdist = z - helix.z0();
+ double safeslope = helix.slope();
+ if (Math.abs(safeslope) < _minslope) safeslope = _minslope * Math.signum(safeslope);
+ return zdist / safeslope;
+ }
+
+ /**
+ * Return a list of x-y path lengths to a cylinder centered on the z axis.
+ * @param helix HelicalTrackFit to be used in calculating the path length
+ * @param r desired radius
+ * @param smax maximum path length to be considered
+ * @param mxint Maximum number of intersections
+ * @return list of path lengths
+ */
+ public static List<Double> PathToCylinder(HelicalTrackFit helix, double r, double smax, int mxint) {
+ // Create a list to hold the path lengths
+ List<Double> pathlist = new ArrayList<Double>();
+ // Retrieve helix dca and RC
+ double dca = helix.dca();
+ double RC = helix.R();
+ // Find cos(dphi) for dphi from the point of closest approach to the first crossing
+ double cdphi = 1 - (r * r - dca * dca) / (2. * RC * (RC - dca));
+ // Make sure we have a physical intersection
+ if (Math.abs(cdphi) < 1.) {
+ // Calculate dphi and the path length to the first crossing
+ double dphi = Math.acos(cdphi);
+ Double s = dphi * Math.abs(RC);
+ // Loop over crossings until we exceed one of the limits
+ while (s < smax && pathlist.size() < mxint) {
+ // Add this odd-numbered crossing to the list
+ pathlist.add(s);
+ // Advance to the next even-numbered crossing
+ s += 2. * (Math.PI - dphi) * Math.abs(RC);
+ // Check to see if we should add it
+ if (s < smax && pathlist.size() < mxint) pathlist.add(s);
+ // Add this even-numbered crossing to the list
+ s += 2. * dphi * Math.abs(RC);
+ }
+ }
+ return pathlist;
+ }
+
+ /**
+ * Return a unit vector giving the track direction at a given point on
+ * the helix.
+ * @param helix HelicalTrackFit to be used in calculating direction
+ * @param s path length to the desired point on helix
+ * @return direction unit vector
+ */
+ public static Hep3Vector Direction(HelicalTrackFit helix, double s) {
+ // Calculate the azimuthal direction
+ double phi = helix.phi0() - s / helix.R();
+ double sth = helix.sth();
+ // Calculate the components of the direction unit vector
+ double ux = Math.cos(phi) * helix.sth();
+ double uy = Math.sin(phi) * helix.sth();
+ double uz = helix.cth();
+ // Return the direction unit vector
+ return new BasicHep3Vector(ux, uy, uz);
+ }
+
+ /**
+ * Return the TrackDirection object for a given point on a helix.
+ * This might seem like an odd place to put this code, but putting
+ * it here allows the cached helix quantities to be used in constructin
+ * the TrackDirection objects.
+ * @param helix HelicalTrackFit to use in constructing the TrackDirection
+ * @param s path length specifying location on helix
+ * @return TrackDirection object for this point on the helix
+ */
+ public static TrackDirection CalculateTrackDirection(HelicalTrackFit helix, double s) {
+ Hep3Vector dir = Direction(helix, s);
+ Matrix deriv = DirectionDerivates(helix, dir, s);
+ return new TrackDirection(dir, deriv);
+ }
+
+ /**
+ * Return the location in space for a particular point on a helix.
+ * @param helix HelicalTrackFit to be used
+ * @param s path length
+ * @return point in space corresponding to the given path length
+ */
+ public static Hep3Vector PointOnHelix(HelicalTrackFit helix, double s) {
+ // Find the azimuthal direction at this path length
+ double RC = helix.R();
+ double phi = helix.phi0() - s / RC;
+ // Calculate the position on the helix at this path length
+ double x = helix.xc() - RC * Math.sin(phi);
+ double y = helix.yc() + RC * Math.cos(phi);
+ double z = helix.z0() + s * helix.slope();
+ // Return the position as a Hep3Vector
+ return new BasicHep3Vector(x, y, z);
+ }
+
+ private static double PathCalc(double xc, double yc, double RC, double x1, double y1, double x2, double y2) {
+ // Find the angle between these points measured wrt the circle center
+ double phi1 = Math.atan2(y1 - yc, x1 - xc);
+ double phi2 = Math.atan2(y2 - yc, x2 - xc);
+ double dphi = phi2 - phi1;
+ // Make sure dphi is in the valid range (-pi, pi)
+ if (dphi > Math.PI) dphi -= 2. * Math.PI;
+ if (dphi < -Math.PI) dphi += 2. * Math.PI;
+ // Return the arc length
+ return -RC * dphi;
+ }
+
+ /**
+ * Return the derivatives of the momentum unit vector with respect to the
+ * helix parameters. The direction derivatives are returned in matrix
+ * form, with the row giving the cartesian component of the momentum
+ * vector and the column giving the helix parameter.
+ * @param helix HelicalTrackFit to be used in calculating derivatives
+ * @param s path length to the desired point on the helix
+ * @return direction derivatives
+ */
+ private static Matrix DirectionDerivates(HelicalTrackFit helix, Hep3Vector u, double s) {
+ // Create the matrix that will hold the derivatives
+ BasicMatrix deriv = new BasicMatrix(3,5);
+ // Retrieve some helix info
+ double cphi0 = Math.cos(helix.phi0());
+ double sphi0 = Math.sin(helix.phi0());
+ double sth = helix.sth();
+ double cth = helix.cth();
+ double dca = helix.dca();
+ double omega = helix.curvature();
+ // Calculate the non-zero derivatives of the direction with respect to the helix parameters
+ deriv.setElement(0, helix.curvatureIndex, (u.x() - cphi0 * sth) / omega); // du_x / domega
+ deriv.setElement(1, helix.curvatureIndex, (u.y() - sphi0 * sth) / omega); // du_y / domega
+ deriv.setElement(0, helix.dcaIndex, -omega * cphi0 * sth); // du_x / ddca
+ deriv.setElement(1, helix.dcaIndex, -omega * sphi0 * sth); // du_y / ddca
+ deriv.setElement(0, helix.phi0Index, -(1 - dca * omega) * sphi0 * sth); // du_x / dphi0
+ deriv.setElement(1, helix.phi0Index, (1 - dca * omega) * cphi0 * sth); // du_y / dphi0
+ deriv.setElement(0, helix.slopeIndex, -u.x() * sth * cth); // du_x / dslope
+ deriv.setElement(1, helix.slopeIndex, -u.y() * sth * cth); // du_y / dslope
+ deriv.setElement(2, helix.slopeIndex, sth*sth*sth); // du_z / dslope
+ // Return the derivative matrix
+ return deriv;
+ }
+/**
+ * return true if the helix is intercepting the bounded cylinder
+ * @param helix
+ * @param r rarius of the cylinder
+ * @param zmin lower bound of the cylinder
+ * @param zmax higher bound of the cylinder
+ * @return
+ */
+ public static boolean isInterceptingBoundedCylinder(HelicalTrackFit helix,double r,double zmin,double zmax){
+ double minpath = PathToZPlane(helix, zmin);//not sure it's very efficient to calculate the maximum path
+ double maxpath = PathToZPlane(helix, zmax);
+ double path = Math.max(minpath, maxpath);
+ List<Double> pathlist = PathToCylinder(helix,r,path,10);
+ for(double s:pathlist){
+ double z = PointOnHelix(helix,s).z();
+ if(z<zmax && z> zmin)
+ return true;
+ }
+ return false;
+ }
+ /**
+ * return true if the helix is intercepting the given disk
+ * @param helix
+ * @param rmin
+ * @param rmax
+ * @param z
+ * @return
+ */
+ public static boolean isInterceptingZDisk(HelicalTrackFit helix,double rmin,double rmax, double z ){
+ double s = PathToZPlane(helix,z);
+ Hep3Vector point = PointOnHelix(helix,s);
+ double x = point.x();
+ double y = point.y();
+ double r = Math.sqrt(x*x+y*y);
+ if(r < rmax && r > rmin ){return true;}
+ return false;
+ }
+
+ /**
+ * return true if the poitn on the helix at coord z is intercepting a given ZPolygon, the
+ * method don't check if the polygone is parallel to a z plane
+ * @param helix
+ * @param z
+ * @param vertices
+ * @return true OR false wether the helix intercept or not the polygone
+ */
+ public static boolean isInterceptingZpolygone(HelicalTrackFit helix,double z,List<Hep3Vector> vertices){
+ //todo: check if all vertices are on same z
+ double epsilon = Math.PI/360.;
+ double s = PathToZPlane(helix,z);
+ if(s<0){return false;}
+ if(s>_maxpath){return false;}
+ Hep3Vector point = PointOnHelix(helix,s);
+ double angle = insidePolygone(point,vertices);
+ if(Math.abs(angle-2*Math.PI)<epsilon)
+ return true;
+ return false;
+ }
+
+
+
+ /**
+ * Check if the given helix is intercepting XY plan
+ * note, the z coordinate of the XY plane should be 0, but due to numerical
+ * approximation, it accept all z value < _epsilon
+ * @param helix the Helix
+ * @param normal a unitary normal vector to the plan
+ * @param orig one point of the plane
+ * @return true if the helix intersect at least once with the plane
+ */
+ public static boolean isInterceptingXYPlane(HelicalTrackFit helix, Hep3Vector normal,Hep3Vector orig){
+ if(normal.z()>_epsilon)
+ throw new UnsupportedOperationException("Not a XY plan !"+normal);
+ double x = -helix.x0()+orig.x();
+ double y = -helix.y0()+orig.y();
+ double z = -helix.z0()+orig.z();
+ double xn= normal.x();
+ double yn= normal.y();
+ double zn= normal.z();
+ double Phip = Math.atan2(yn, xn);
+ double dist = (xn*x+yn*y+zn*z);
+ double verif = Math.sin(Phip-helix.phi0())+helix.curvature()*dist;
+ if(normal.magnitude() < 0.99)
+ throw new UnsupportedOperationException("normal vector not unitary :"+normal.magnitude());
+ if(Math.abs(verif)>1)
+ return false;
+ else
+ return true ;
+ }
+
+ /**
+ * Check if the given helix is intercepting an XY plane bouded by a
+ * list of nodes
+ * @param helix
+ * @param normal normal vector to the plane
+ * @param nodes
+ * @return true if th e helix intercept the bounded plane
+ **/
+ public static boolean isInterceptingBoundedXYPlane(HelicalTrackFit helix,Hep3Vector normal,List<Hep3Vector> nodes){
+ Hep3Vector orig = nodes.get(0);
+ if(!isInterceptingXYPlane(helix,normal,orig)){return false;}
+ double s = PathToXYPlan(helix, normal, orig);
+ if(s<0) {return false;}
+ if(s>_maxpath){return false;}
+ Hep3Vector point = PointOnHelix(helix, s);
+ if(Math.abs(insidePolygone(point, nodes)-2*Math.PI)<_epsilon)
+ return true;
+ return false;
+ }
+ /**
+ * return one path length (projection on the XY plane) to an XY plane
+ * the methode is for now only used to check IF an helix is intercepting an
+ * XY plane
+ * note, the z coordinate of the XY plane should be 0, but due to numerical
+ * approximation, it accept all z value < _epsilon
+ * @param helix
+ * @param normal a UNITARY vector NORMAL to the plan
+ * @param origin one point of the plane
+ * @return path length
+ */
+ public static double PathToXYPlan(HelicalTrackFit helix,Hep3Vector normal,Hep3Vector origin){
+ if(normal.z()>_epsilon)
+ throw new UnsupportedOperationException("Not a XY plan !"+normal);
+ double x = -helix.x0()+origin.x();
+ double y = -helix.y0()+origin.y();
+ double z = -helix.z0()+origin.z();
+ double xn= normal.x();
+ double yn= normal.y();
+ double zn= normal.z();
+ double phip = Math.atan2(yn, xn);
+ double dist = (xn*x+yn*y+zn*z);
+ double sinPlan = Math.sin(phip-helix.phi0())+helix.curvature()*dist;
+ double Cs = (Math.asin(sinPlan)-phip+helix.phi0());
+ double s=dist/(Math.sqrt(xn*xn+yn*yn)*(2./Cs)*Math.sin(Cs/2.)*Math.cos(Cs/2+phip-helix.phi0())+zn*helix.slope());
+ return s;
+ }
+
+
+ /**
+ * adapted from http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/
+ * find if a point is inside a given polygon (!!! not polyhedron !!!) in 3Dspace
+ * if the point is in the polygone the returned value is 2*PI
+ * if the point is outside of the polygone the returned value is 0
+ * if the point is not of the plan of the polygone the value is between 0 and 2PI
+ * The polygone HAVE to Be convex for the algorythme to work
+ *
+ * This fonction sum all the angle made by q and 2 consecutives point of the polygone
+ *
+ * @param q the point to determine if it is or not in the polygon
+ * @param p List of edges of the polygone
+ * @return 2PI if inside, 0 if outside,intermediate values if not on the plane
+ */
+ private static double insidePolygone(Hep3Vector q,List<Hep3Vector> p)
+ {
+ int i;
+ double Epsilon = _epsilon;
+ double m1,m2;
+ double anglesum=0,costheta;
+ Hep3Vector p1;
+ Hep3Vector p2;
+ int n = p.size();
+ for (i=0;i<n;i++) {
+ p1 = p.get(i);
+ p2 = p.get((i+1)%n);
+ double x1 = (p1.x()-q.x());
+ double y1 = (p1.y()-q.y());
+ double z1 = (p1.z()-q.z());
+ double x2 = (p2.x()-q.x());
+ double y2 = (p2.y()-q.y());
+ double z2 = (p2.z()-q.z());
+ m1 = Math.sqrt(x1*x1+y1*y1+z1*z1);
+ m2 = Math.sqrt(x2*x2+y2*y2+z2*z2);
+ if (m1*m2 <= Epsilon){
+ return 2*Math.PI; /* We are on a node, consider this inside */
+ }
+ else{
+ costheta = (x1*x2 + y1*y2 + z1*z2)/(m1*m2);
+ anglesum += Math.acos(costheta);
+ }
+ }
+ return anglesum;
+ }
+
+
+ private static double RC(CircleFit cfit) {
+ return 1.0 / cfit.curvature();
+ }
+
+ private static double xc(CircleFit cfit) {
+ return cfit.xref() + (RC(cfit) - cfit.dca()) * Math.sin(cfit.phi());
+ }
+
+ private static double yc(CircleFit cfit) {
+ return cfit.yref() - (RC(cfit) - cfit.dca()) * Math.cos(cfit.phi());
+ }
+
+ private static double x0(CircleFit cfit) {
+ return cfit.xref() - cfit.dca() * Math.sin(cfit.phi());
+ }
+
+ private static double y0(CircleFit cfit) {
+ return cfit.yref() + cfit.dca() * Math.cos(cfit.phi());
+ }
+
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -N JetDriverJetCount.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ JetDriverJetCount.java 19 Jun 2009 19:11:19 -0000 1.1
@@ -0,0 +1,151 @@
+package org.lcsim.contrib.Mbussonn;
+
+import hep.aida.IHistogramFactory;
+import hep.aida.IProfile1D;
+import hep.physics.jet.JadeEJetFinder;
+import hep.physics.jet.JetFinder;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.HepLorentzVector;
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.HashMap;
+import java.util.LinkedList;
+import java.util.List;
+import java.util.Map;
+import org.lcsim.digisim.MyLCRelation;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.Vertex;
+import org.lcsim.event.base.MCReconstructedParticle;
+import org.lcsim.event.util.JetDriver;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ * A simple driver which can be used to find jets from ReconstructedParticles.
+ * The resuslting jets are stored in a new collection of ReconstructedParticles.
+ * @author tonyj
+ * @version $Id: JetDriverJetCount.java,v 1.1 2009/06/19 19:11:19 mbussonn Exp $
+ */
+public class JetDriverJetCount extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ private IProfile1D g1,g2;
+ private static final Hep3Vector IP = new BasicHep3Vector(0,0,0);
+ private static final String defaultOutputCollectionName = "Jets";
+ private String inputCollectionName;
+ private String outputCollectionName = defaultOutputCollectionName;
+ private JetFinder finder;
+ List<LCRelation> rc2mc = new ArrayList<LCRelation>();
+ List<ReconstructedParticle> collectionOfRCP;
+
+ /** Creates a new instance of JetFinder with the default properties */
+ public JetDriverJetCount()
+ {
+ IHistogramFactory hf = aida.histogramFactory();
+ g1 = hf.createProfile1D("njet vs log2ref",16,0.,16.);
+ g2 = hf.createProfile1D("njet vs ref",110,0.,.5);
+
+ }
+
+
+ public String getInputCollectionName()
+ {
+ return inputCollectionName;
+ }
+
+ /**
+ * The name of the input collection to use. If not set (or set to <code>null</code> uses
+ * the first collection of ReconstructedParticles found in the event.
+ * @param inputCollectionName
+ */
+
+ public void setInputCollectionName(String inputCollectionName)
+ {
+ this.inputCollectionName = inputCollectionName;
+ }
+
+ public String getOutputCollectionName()
+ {
+ return outputCollectionName;
+ }
+
+ /**
+ * The name of the output collection added to the event. If not set, or set to <code>null</code>,
+ * default to "Jets"
+ * @param outputCollectionName
+ */
+
+ public void setOutputCollectionName(String outputCollectionName)
+ {
+ this.outputCollectionName = outputCollectionName == null ? defaultOutputCollectionName : outputCollectionName;
+ }
+
+ public JetFinder getFinder()
+ {
+ return finder;
+ }
+
+ /**
+ * Set the jet finding algorithm to use
+ * @param finder
+ */
+ @Override
+ protected void process(EventHeader event)
+ {
+ super.process(event);
+ // Find the input reconstructed Particles
+ List<ReconstructedParticle> input = null;
+ collectionOfRCP = new ArrayList<ReconstructedParticle>();
+ if (inputCollectionName ==null)
+ {
+ List<List<ReconstructedParticle>> listOfLists = event.get(ReconstructedParticle.class);
+ if (listOfLists.isEmpty()) {
+ input = new LinkedList<ReconstructedParticle>();
+ //List<LCRelation> mcrelations = event.get(LCRelation.class, "HelicalTrackMCRelations");
+
+ for(MCParticle mcp : event.getMCParticles()){
+ if(mcp.getGeneratorStatus()==MCParticle.FINAL_STATE){
+ MCReconstructedParticle rcp = new MCReconstructedParticle(mcp); // reconstruct particle from mcp
+ input.add(rcp); // for further treatment in this driver
+ collectionOfRCP.add(rcp); // list of rcp written line 211
+ }
+ }
+ }
+ else{
+ input = listOfLists.get(0);
+ }
+ }
+ else
+ {
+ input = event.get(ReconstructedParticle.class,inputCollectionName);
+ }
+ Map<HepLorentzVector, ReconstructedParticle> map = new HashMap<HepLorentzVector, ReconstructedParticle>();
+ for (ReconstructedParticle p : input)
+ {
+ map.put(p.asFourVector(),p);
+ }
+ // Pass the data to the Jet finder
+ //finder.setEvent(map.keySet());
+
+
+ int i;
+ for(i=0;i<16;i++)
+ {
+ double res = 1./Math.pow(2,i);
+ finder = new JadeEJetFinder(res);
+ finder.setEvent(map.keySet());
+
+ int njet = finder.njets();
+ System.out.println("njets"+njet+" i:"+i+" res:"+res);
+ aida.cloud2D("number of jet vs res ").fill(res,njet);
+ aida.cloud2D("number of jet vs -log2(res) ").fill(i,njet);
+ g1.fill(i,njet);
+ g2.fill(res,njet);
+ }
+ }
+
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -u -r1.2 -r1.3
--- JetDriverExtended.java 8 Jun 2009 18:44:44 -0000 1.2
+++ JetDriverExtended.java 19 Jun 2009 19:11:19 -0000 1.3
@@ -26,7 +26,7 @@
* A simple driver which can be used to find jets from ReconstructedParticles.
* The resuslting jets are stored in a new collection of ReconstructedParticles.
* @author tonyj
- * @version $Id: JetDriverExtended.java,v 1.2 2009/06/08 18:44:44 mbussonn Exp $
+ * @version $Id: JetDriverExtended.java,v 1.3 2009/06/19 19:11:19 mbussonn Exp $
*/
public class JetDriverExtended extends Driver
{
@@ -92,15 +92,13 @@
protected JetFinder defaultJetFinder()
{
- return new JadeEJetFinder(0.045);
+ return new JadeEJetFinder(0.002);
}
@Override
protected void process(EventHeader event)
{
- super.process(event);
-
- boolean hist = getHistogramLevel() > 0;
+ super.process(event);
double totalEnergy=0;
// Find the input reconstructed Particles
List<ReconstructedParticle> input = null;
@@ -119,9 +117,8 @@
MyLCRelation lcr = new MyLCRelation(rcp, mcp); // new relation ship
rc2mc.add(lcr); // array of relationship written on this file line 213
collectionOfRCP.add(rcp); // list of rcp written line 211
- System.out.println("lcr to string mcp :"+lcr.getTo().toString()+" rcp "+lcr.getFrom().getClass());
+ //System.out.println("lcr to string mcp :"+lcr.getTo().toString()+" rcp "+lcr.getFrom().getClass());
totalEnergy += mcp.getEnergy();
-
}
}
}
@@ -169,7 +166,7 @@
aida.cloud1D("One particle/jet |H4|").fill(finder.jet(i).magnitude());
aida.cloud1D("One particle/jet pr").fill(pr);
aida.cloud1D("One particle/jet pz").fill(pz);
- aida.cloud2D("One particle/jet direction").fill(pz,pr);
+ //aida.cloud2D("One particle/jet direction").fill(pz,pr);
}else{
@@ -179,7 +176,7 @@
aida.cloud1D("Jet particle/jet |H4|").fill(finder.jet(i).magnitude());
aida.cloud1D("Jet particle/jet pr").fill(pr);
aida.cloud1D("Jet particle/jet pz").fill(pz);
- aida.cloud2D("Jet particle/jet direction").fill(pz,pr);
+ //aida.cloud2D("Jet particle/jet direction").fill(pz,pr);
aida.cloud2D("jet direction ").fill(hlv.v3().x(),hlv.v3().y());
aida.cloud1D("JetDriver/particlesPerJet").fill(finder.nParticlesPerJet(i));