5 added files
lcd/hep/lcd/recon/cluster/MIPMatching
diff -N AbstractChi2.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ AbstractChi2.java 10 Aug 2005 17:38:17 -0000 1.1
@@ -0,0 +1,41 @@
+package hep.lcd.recon.cluster.MIPMatching;
+
+import java.io.*;
+import java.util.*;
+
+/**
+ AbstractChi2 for MIP Matching. The errors on the track parameters
+ are derived from studies of single pion samples. These can be
+ overwritten in user defined Chi2-Classes.
+
+ @version 1.0
+ @author [log in to unmask]
+
+ Modification Log:
+ -- 08/08/2005 Version 1.0
+ Calculation of chi2 separated from MIPFit to provide user with
+ possibility to define individual chi2-criteria.
+
+*/
+
+public abstract class AbstractChi2{
+ public AbstractChi2(){}
+
+ public void registerTrackParameters(TrackParameters tpin1, TrackParameters tpin2){
+ tp1 = tpin1;
+ tp2 = tpin2;
+ }
+
+ public abstract double getChi2();
+
+ protected TrackParameters tp1 = null;
+ protected TrackParameters tp2 = null;
+//
+// Typical uncertainty (sqrt(var)) for track parameters from singe pions
+//
+ protected double ed0 = 2.6;
+ protected double ephi0 = 0.06;
+ protected double ekappa = 0.0005;
+ protected double etanl = 0.09;
+ protected double ez0 = 1.;
+}
lcd/hep/lcd/recon/cluster/MIPMatching
diff -N MIPFit.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MIPFit.java 10 Aug 2005 17:38:17 -0000 1.1
@@ -0,0 +1,397 @@
+package hep.lcd.recon.cluster.MIPMatching;
+
+import java.io.*;
+import java.util.*;
+
+import hep.lcd.event.*;
+import hep.lcd.geometry.CalorimeterCell;
+import hep.lcd.util.driver.AbstractProcessor;
+
+import hep.physics.*;
+import hep.lcd.recon.cluster.MIP.*;
+/**
+ Driver for MIP Track Matching.
+ <p>
+ Example of Use:
+ <p><blockquote><pre>
+ MIPFit mipFit = new MIPFit();
+ AbstractChi2 chi2 = new ThreeDChi2();
+
+ mipFit.setTag("MIPCluster EMCal");
+ mipFit.registerChi2(chi2);
+ mipFit.setChi2Max(40);
+
+ add(mipFit);
+ </pre></blockquote><p>
+ The MIPClusterList is read from <tt>LCDEvent</tt>. The name under
+ which it is stored can be defined by the user via
+ <tt>mipFit.setTag(String)</tt>. The same tag is used to write out the
+ new <tt>MIPClusterList</tt> containing merged objects.
+
+ The Helix parameters for each MIPCluster is calculated using the
+ IP, the first and the last hit of the corresponding MIPCluster in
+ the EMCal.
+
+ The Chi2 criteria for the matching of MIPClusters can be defined by
+ the user. The corresponding class has to be derived from
+ <tt>AbstractChi2.java</tt> and has to provide the method
+ <tt>getChi2()</tt>. As a default, the class <tt>ThreeDChi2</tt> is
+ provided, assuming typical errors for the helix parameters
+ estimated from studies based on single pion samples. The default
+ <tt>chi2</tt> for the matching of MIPClusters is <tt>40</tt>.
+
+ @author Wolfgang F. Mader ([log in to unmask])
+ @version 1.0
+
+ Modification Log:
+ -- 07/22/2004 Version 1.0
+ -- 08/08/2005 Version 1.1
+ Addapted for use with general MIP Tracking Package
+*/
+
+public class MIPFit extends AbstractProcessor{
+//
+// Field Summary
+//
+ private CalorimeterCell thisCell = null;
+ private boolean lDebug = true;
+ private Vector usedMIPClusters;
+ private Vector vMIPClusters = null;
+ private String objectTag = new String("MIPCluster");
+
+ private AbstractChi2 chi2 = null;
+ private double chi2Max = 40;
+
+//
+// Constructor Summary
+//
+ public MIPFit(){}
+
+ public void init(){}
+ public void start(){}
+ public void stop(){}
+//
+// Method Summary
+//
+ public void setChi2Max(double d){
+ chi2Max = d;
+ }
+
+ public void registerChi2(AbstractChi2 c){
+ chi2 = c;
+ }
+
+ public void setTag(String s){
+ objectTag = s;
+ }
+
+ public void setlDebug(boolean b){
+ lDebug = b;
+ }
+//
+// Main Process Loop
+//
+ public void process(LCDEvent event){
+
+ thisEvent = event;
+ thisCell = thisEvent.getCalorimeterCell();
+//
+// No Assertion Error needed -- thrown by LCDEvent
+//
+ MIPClusterList lMIPClusters = (MIPClusterList) thisEvent.get(objectTag);
+
+ if ( lDebug )
+ System.out.println("Number of MIP Objects found "+lMIPClusters.getNClusters());
+ if ( lMIPClusters.getNClusters() < 2 ) return;
+
+ initialize(lMIPClusters.getClusters());
+
+ Vector vMIPObjects = new Vector();
+ Vector vMIPPairs = new Vector();
+//
+// Calculate Chi2 for all MIP combinations
+//
+ for ( int i=0; i<vMIPClusters.size()-1; i++) {
+ for ( int j=i+1; j<vMIPClusters.size(); j++) {
+
+ System.out.println("i "+i+" j "+j);
+ MIPCluster mip1 = (MIPCluster) vMIPClusters.elementAt(i);
+ MIPCluster mip2 = (MIPCluster) vMIPClusters.elementAt(j);
+
+ vMIPPairs.add(pairObjects(mip1,mip2));
+ }
+ }
+
+ if ( lDebug ) {
+ Enumeration e = vMIPPairs.elements();
+ while ( e.hasMoreElements() ){
+ MIPPair mipPair = (MIPPair) e.nextElement();
+ System.out.println("1: "+mipPair.getFirstMIP()+
+ " 2: "+mipPair.getSecondMIP()+
+ " chi2 "+mipPair.getChi2());
+ }
+ }
+//
+// Sort MIPs according to chi2
+//
+ Collections.sort(vMIPPairs,new MIPChi2Sort());
+
+ if ( lDebug ) {
+ Enumeration e = vMIPPairs.elements();
+ while ( e.hasMoreElements() ){
+ MIPPair mipPair = (MIPPair) e.nextElement();
+ System.out.println("--1: "+mipPair.getFirstMIP()+" "+mipPair.getFirstMIP().getNHits()+
+ " 2: "+mipPair.getSecondMIP()+" "+mipPair.getSecondMIP().getNHits()+
+ " chi2 "+mipPair.getChi2());
+ }
+ }
+
+ Enumeration e = vMIPPairs.elements();
+ while ( e.hasMoreElements() ){
+ MIPPair mipPair = (MIPPair) e.nextElement();
+
+ MIPCluster mip1 = mipPair.getFirstMIP();
+ MIPCluster mip2 = mipPair.getSecondMIP();
+
+ int index1 = vMIPClusters.indexOf(mip1);
+ int index2 = vMIPClusters.indexOf(mip2);
+
+ boolean b1 = ((Boolean) usedMIPClusters.elementAt(index1)).booleanValue();
+ boolean b2 = ((Boolean) usedMIPClusters.elementAt(index2)).booleanValue();
+
+//
+// Pick MIPCluseter combinations as long as chi2 < chi2max...
+//
+ if ( !(b1 || b2) && mipPair.getChi2()<chi2Max ){
+//
+// ... and merge the corresponding MIPClusters
+//
+ Enumeration eHits = mip2.getHits();
+ while ( eHits.hasMoreElements() ){
+ CalorimeterHit hit = (CalorimeterHit) eHits.nextElement();
+ mip1.addHit(hit);
+ }
+ vMIPObjects.add(mip1);
+//
+// set 'used' flag
+//
+ usedMIPClusters.removeElementAt(index1);
+ usedMIPClusters.insertElementAt(new Boolean(true),index1);
+ usedMIPClusters.removeElementAt(index2);
+ usedMIPClusters.insertElementAt(new Boolean(true),index2);
+ } else {
+//
+// Fill remaining MIPClusters to output vector and set 'used' flag
+//
+ if ( !b1 ){
+ vMIPObjects.add(mip1);
+ usedMIPClusters.removeElementAt(index1);
+ usedMIPClusters.insertElementAt(new Boolean(true),index1);
+ }
+ if ( !b2 ){
+ vMIPObjects.add(mip2);
+ usedMIPClusters.removeElementAt(index2);
+ usedMIPClusters.insertElementAt(new Boolean(true),index2);
+ }
+ }
+ }
+
+ if ( lDebug ) {
+ Enumeration en = vMIPObjects.elements();
+ while ( en.hasMoreElements() ) {
+ MIPCluster cluster = (MIPCluster) en.nextElement();
+ System.out.println("clusters "+cluster+" "+cluster.getNHits());
+ }
+ }
+//
+// Write information out the LCDEvent
+//
+ thisEvent.put(objectTag,new MIPClusterList(vMIPObjects));
+ }
+
+ private MIPPair pairObjects(MIPCluster mip1, MIPCluster mip2){
+
+//
+// Initialize Helix calculation
+//
+ MIPHelix h1 = new MIPHelix(mip1,thisCell);
+ MIPHelix h2 = new MIPHelix(mip2,thisCell);
+
+//
+// Get Helix Parameters
+//
+ if ( lDebug ){
+ System.out.println("MIP1: First Hit -- "+getFirstHitIndex(mip1));
+ System.out.println(" Last Hit -- "+getLastHitIndex(mip1));
+ }
+ h1.getEstimate(-1,getFirstHitIndex(mip1),getLastHitIndex(mip1));
+
+ if ( lDebug ){
+ System.out.println("MIP2: First Hit -- "+getFirstHitIndex(mip2));
+ System.out.println(" Last Hit -- "+getLastHitIndex(mip2));
+ }
+ h2.getEstimate(-1,getFirstHitIndex(mip2),getLastHitIndex(mip2));
+
+//
+// Store Information in TrackParameters Object
+//
+ TrackParameters tp1 = h1.getTrackParameters();
+ if ( lDebug ) printTrackParameters(tp1);
+
+ TrackParameters tp2 = h2.getTrackParameters();
+ if ( lDebug ) printTrackParameters(tp2);
+
+//
+// Tell chi2 Caluclator about track parameters
+//
+ chi2.registerTrackParameters(tp1, tp2);
+//
+// Store information in internal class MIPPair for further processing
+//
+ MIPPair mipPair = new MIPPair(mip1,mip2,tp1,tp2,chi2.getChi2());
+ return mipPair;
+ }
+
+//
+// Get index of innermost layer
+//
+ private int getFirstHitIndex(MIPCluster mip){
+ Enumeration e = mip.getHits();
+
+ int firstLayer = 100;
+ int firstHitIndex = -1;
+
+ int i=0;
+ while ( e.hasMoreElements() ){
+ CalorimeterHit hit = (CalorimeterHit) e.nextElement();
+ thisCell.setTowerID(hit.getTowerID());
+
+ int thisLayer = thisCell.getLayer();
+ if ( thisLayer < firstLayer ) {
+ firstLayer = thisLayer;
+ firstHitIndex = i;
+ }
+ i++;
+ }
+ return firstHitIndex;
+ }
+
+//
+// Get index of outermost layer
+//
+ private int getLastHitIndex(MIPCluster mip){
+ Enumeration e = mip.getHits();
+
+ int lastLayer = 0;
+ int lastHitIndex = 0;
+
+ int i=0;
+ while ( e.hasMoreElements() ){
+ CalorimeterHit hit = (CalorimeterHit) e.nextElement();
+ thisCell.setTowerID(hit.getTowerID());
+
+ int thisLayer = thisCell.getLayer();
+ if ( thisLayer > lastLayer ){
+ lastLayer = thisLayer;
+ lastHitIndex = i;
+ }
+ i++;
+ }
+ return lastHitIndex;
+ }
+//
+// Store MIP information from LCDEvent in temporary vector for further processing
+//
+ private void initialize(Enumeration e){
+ int i=0;
+ vMIPClusters = new Vector();
+ usedMIPClusters = new Vector();
+ while ( e.hasMoreElements() ){
+ System.out.println("i"+i++);
+ vMIPClusters.add((MIPCluster) e.nextElement());
+ usedMIPClusters.add(new Boolean(false));
+ }
+ }
+//
+// Print out track parameters ( for debugging )
+//
+ private void printTrackParameters(TrackParameters tp){
+ System.out.println("kappa "+tp.kappa());
+ System.out.println("d0 "+tp.d0());
+ System.out.println("phi0 "+tp.phi0());
+ System.out.println("z0 "+tp.z0());
+ System.out.println("tanl "+tp.tanl());
+ System.out.println();
+ }
+//
+// Field Summary
+//
+ private LCDEvent thisEvent = null;
+//
+// Track Parameters
+//
+}
+//
+// Internal convenience class MIPPair for information processing
+//
+class MIPPair{
+ MIPPair(MIPCluster mip1,
+ MIPCluster mip2,
+ TrackParameters tp1,
+ TrackParameters tp2,
+ double in){
+ m1 = mip1;
+ m2 = mip2;
+
+ trackParameters1 = tp1;
+ trackParameters2 = tp2;
+
+ chi2 = in;
+ }
+
+// public MIPCluster[] getMIPCluster(){
+// MIPCluster
+// }
+
+ public double getChi2(){
+ return chi2;
+ }
+
+ public MIPCluster getFirstMIP(){
+ return m1;
+ }
+ public MIPCluster getSecondMIP(){
+ return m2;
+ }
+
+ MIPCluster m1;
+ MIPCluster m2;
+ double chi2;
+ TrackParameters trackParameters1;
+ TrackParameters trackParameters2;
+}
+//
+// Comparator class for Chi2-Sorting
+//
+class MIPChi2Sort implements Comparator{
+
+ /** The compare function used for sorting.
+ * Comparison is done on cluster energy.
+ * @param obj1 MIPPair1
+ * @param obj2 MIPPair2
+ * @return
+ * <ol>
+ * <li> -1 if MIPPair1 > MIPPair2
+ * <li> 0 if MIPPair1 = MIPPair2
+ * <li> 1 if MIPPair1 < MIPPair2
+ * </ol>
+ */
+ public int compare(Object obj1, Object obj2)
+ {
+ if(obj1==obj2) return 0;
+ MIPPair mip1 = (MIPPair) obj1;
+ MIPPair mip2 = (MIPPair) obj2;
+ if(mip1.getChi2()-mip2.getChi2()<0.) return -1;
+ return 1;
+ }
+}
lcd/hep/lcd/recon/cluster/MIPMatching
diff -N MIPHelix.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MIPHelix.java 10 Aug 2005 17:38:18 -0000 1.1
@@ -0,0 +1,397 @@
+package hep.lcd.recon.cluster.MIPMatching;
+import hep.lcd.recon.cluster.MIP.*;
+
+import java.util.Enumeration;
+
+import hep.lcd.geometry.CalorimeterCell;
+import hep.lcd.event.CalorimeterHit;
+
+import hep.lcd.event.LCDEvent;
+import hep.lcd.event.TrackerHits;
+import hep.lcd.event.TrackerHit;
+/**
+ MIPHelix. Calculates Estimate of Helix Parameters.
+
+ @author Wolfgang F. Mader ([log in to unmask])
+ @version 1.0
+
+ Modification Log:
+ -- 07/22/2004 Version 1.0
+*/
+
+public class MIPHelix{
+//
+// Constructor Summary
+//
+ MIPHelix (MIPCluster mipObject, CalorimeterCell cell){
+
+ thisCell = cell;
+ thisMIPObject = mipObject;
+ registerHitSamples();
+
+ tP = new TrackParameters();
+ }
+
+//
+// Method Summary
+//
+ private void registerHitSamples(){
+
+ Enumeration e = thisMIPObject.getHits();
+
+ int i=0;
+ System.out.println("-->"+i);
+ while ( e.hasMoreElements() ){
+ CalorimeterHit hit = (CalorimeterHit) e.nextElement();
+ thisCell.setTowerID(hit.getTowerID());
+
+ double[] tmp = thisCell.getPosition();
+ for ( int j=0; j<3; j++) Hits[i][j] = tmp[j];
+
+ errorHits[i][0] = 0.5/Math.sqrt(12);
+ errorHits[i][1] = 0.5/Math.sqrt(12);
+ errorHits[i][2] = 0.;
+
+ i++;
+ }
+ }
+
+ public void getEstimate(int i, int j, int k){
+ int lastHit = numberOfHits-1;
+ if ( i == -1 ){
+ X = new double[] {0.0001, Hits[j][0], Hits[k][0]};
+ Y = new double[] {0.0001, Hits[j][1], Hits[k][1]};
+ Z = new double[] {0.0001, Hits[j][2], Hits[k][2]};
+ } else {
+ X = new double[] {Hits[i][0], Hits[j][0], Hits[k][0]};
+ Y = new double[] {Hits[i][1], Hits[j][1], Hits[k][1]};
+ Z = new double[] {Hits[i][2], Hits[j][2], Hits[k][2]};
+ }
+ tP.updateParameters(kappa(), d0(), phi0(), z0(), lambda());
+ }
+
+ public TrackParameters getTrackParameters(){
+ return tP;
+ }
+
+ private double arcLength(double X, double Y){
+ double length = -1;
+
+ double AKAD = tP.kappa();
+ double ALPHAD = tP.alpha();
+ double BETAD = tP.beta();
+ double GAMMAD = tP.gamma();
+ double XID = Math.sqrt(1.0+4*AKAD*GAMMAD);
+ double R2KAP = 2.0 * AKAD;
+ double XINV = 1./XID;
+ double AR2KAP = Math.abs(R2KAP);
+
+ double DELH1 = AKAD*(Math.pow(X,2) + Math.pow(Y,2)) + ALPHAD*X + BETAD*Y + GAMMAD;
+ double DELH = DELH1*(1. - AKAD*DELH1);
+//
+// use an exact calculation if we need to...
+//
+ double AKADEL = 4.*AKAD *DELH1;
+ if ( Math.abs(DELH*AKAD) > 0.01 && AKADEL > -1.) {
+ DELH = 0.5/AKAD*( Math.sqrt(1.+ AKADEL) - 1.);
+ }
+ double XON = (X - DELH*ALPHAD)/(1.+R2KAP*DELH);
+ double YON = (Y - DELH*BETAD )/(1.+R2KAP*DELH);
+
+ double CROSS = ALPHAD*YON - BETAD*XON;
+ double A1 = AR2KAP*CROSS*XINV;
+ double A2 = R2KAP*(ALPHAD*XON + BETAD*YON)*XINV + XID;
+
+ if (A1 >= 0.0 && A2 > 0.0 && A1 < 0.3){
+ double ARG2 = Math.pow(A1,2);
+ length = CROSS*(1. + ARG2*(1./6. + ARG2*(3./40.)))*XINV;
+ } else {
+ length = Math.atan2(A1,A2);
+ if (length < 0.) length += 2.*Math.PI;
+ length = length/AR2KAP;
+ }
+ return length;
+ }
+
+ /*
+Formulas for fitting the stuff in z-direction
+ */
+ private double Xij(int i, int j){
+ return ((double)
+ X[i-1]-X[j-1]
+ );
+ }
+
+ private double Yij(int i, int j){
+ return ((double)
+ Y[i-1]-Y[j-1]
+ );
+ }
+
+ private double Delta(){
+ return ((double)
+ Deltaij(1,2)+Deltaij(2,3)+Deltaij(3,1)
+ );
+ }
+
+ private double Deltaij(int i, int j){
+ return ((double)
+ X[i-1]*Y[j-1] - X[j-1]*Y[i-1]
+ );
+ }
+
+ private double R(int i){
+ return ((double) Math.sqrt(X[i-1]*X[i-1] + Y[i-1]*Y[i-1]) );
+ }
+
+ private double Rij(int i, int j){
+ return ((double)
+ R(i) - R(j)
+ );
+ }
+
+ private double kappa(){
+ return ((double)
+ Delta()/(Rij(1,2)*Rij(2,3)*Rij(3,1))
+ );
+ }
+
+ private double prime(String ab){
+ double XY31 = 0;
+ double XY21 = 0;
+
+ if ( ab.equals("alpha") ){
+ XY31 = Yij(3,1);
+ XY21 = -Yij(2,1);
+ } else if ( ab.equals("beta") ) {
+ XY31 = -Xij(3,1);
+ XY21 = Xij(2,1);
+ }
+ return ((double) -kappa()/Delta()
+ *( XY31 * Rij(2,1)*Rij(2,1)
+ +XY21 * Rij(3,1)*Rij(3,1) )
+ );
+ }
+
+ private double gamma(){
+ return ((double)
+ kappa()*R(1)*R(1)
+ - prime("alpha")*X[0]
+ - prime("beta")*Y[0]
+ );
+ }
+
+ private double alpha(){
+ double dkappa = kappa();
+ double dgamma = gamma();
+ double dR1 = R(1);
+ double dR2 = R(2);
+
+ return ((double)
+ -1/Deltaij(1,2) * ( (dkappa*dR1*dR1+dgamma)*Y[1]
+ -(dkappa*dR2*dR2+dgamma)*Y[0] )
+ );
+ }
+
+ private double beta(){
+ double dkappa = kappa();
+ double dgamma = gamma();
+ double dR1 = R(1);
+ double dR2 = R(2);
+
+ return ((double)
+ -1/Deltaij(1,2) * ( -(dkappa*dR1*dR1+dgamma)*X[1]
+ +(dkappa*dR2*dR2+dgamma)*X[0] )
+ );
+ }
+
+ private double xi(){
+ double dalpha = alpha();
+ double dbeta = beta();
+ return ((double)
+ Math.sqrt(dalpha*dalpha + dbeta*dbeta)
+ );
+ }
+
+ private double xikappagamma(){
+ return ((double)
+ Math.sqrt(1+4*kappa()*gamma())
+ );
+ }
+
+ private double d0(){
+ return ((double)
+ (xi()-1.0)/(2*kappa())
+ );
+ }
+
+ private double d0gamma(){
+ double dkappa = kappa();
+ return ((double)
+ Math.sqrt( gamma()/dkappa + 1/(4*dkappa*dkappa) )
+ - 1/(2*dkappa)
+ );
+ }
+
+ private double d0xi(){
+ return ((double)
+ 2*gamma()/(xi()+1)
+ );
+ }
+
+ private double phi0(){
+ return ((double)
+ Math.atan2(alpha(),-beta())
+ );
+ }
+
+ /*
+Formulas for fitting the stuff in z-direction
+ */
+
+ private double d(double xIn, double yIn){
+ return ((double)
+ (Math.sqrt( Math.max(1+4*kappa()*delta(xIn, yIn),0) ) - 1)
+ /( 2*kappa() )
+ );
+ }
+
+ private double dddkappa(double xIn, double yIn){
+ double dkappa = kappa();
+ double dalpha = alpha();
+ double dbeta = beta();
+ double dgamma = gamma();
+ double cX = CX(xIn, yIn);
+
+ double dtemp = 1./(2.*dkappa*dkappa) *
+ (
+ 1+1./(2.*cX) * ( -1/dkappa - 2.*(dalpha*xIn+dbeta*yIn+dgamma))
+ );
+
+ return((double) dtemp);
+ }
+
+ private double dddabg(double xIn, double yIn, double xyIn){
+ double cX = CX(xIn, yIn);
+ double dkappa = kappa();
+ double dtemp = 1./(2.*cX) * (xyIn/dkappa);
+
+ return((double) dtemp);
+ }
+
+ private double CX(double xIn, double yIn){
+ double dkappa = kappa();
+ double dalpha = alpha();
+ double dbeta = beta();
+ double dgamma = gamma();
+
+ double dtemp = 1./(4.*dkappa*dkappa) +
+ 1./dkappa * ( dalpha*xIn + dbeta*yIn + dgamma ) +
+ Math.sqrt(xIn*xIn + yIn*yIn);
+
+ return ((double) dtemp);
+ }
+
+ private double delta(double xIn, double yIn){
+ return ((double)
+ kappa()*(xIn*xIn + yIn*yIn)
+ + alpha()*xIn + beta()*yIn + gamma()
+ );
+ }
+
+ private double xOn(double xIn, double yIn){
+ double dd = d(xIn, yIn);
+
+ double Nominator = xIn - dd*alpha();
+ double deNominator = 1 + 2*kappa()*dd;
+ if ( deNominator < 0. ) deNominator = Math.max(-400,deNominator);
+ if ( deNominator > 0. ) deNominator = Math.min( 400,deNominator);
+ if ( deNominator == 0. ) deNominator = 0.00001;
+
+ return ((double)
+ Nominator / deNominator
+ );
+ }
+
+ private double yOn(double xIn, double yIn){
+ double dd = d(xIn, yIn);
+
+ double Nominator = yIn - dd*beta();
+ double deNominator = 1 + 2*kappa()*dd;
+
+ if ( deNominator < 0. ) deNominator = Math.max(-400,deNominator);
+ if ( deNominator > 0. ) deNominator = Math.min( 400,deNominator);
+ if ( deNominator == 0. ) deNominator = 0.00001;
+
+ return ((double)
+ Nominator / deNominator
+ );
+ }
+
+ private double cross(double xIn, double yIn){
+ return ((double)
+ alpha()*yOn(xIn,yIn) - beta()*xOn(xIn,yIn)
+ );
+ }
+
+ private double aOne(double xIn, double yIn){
+ return ((double)
+ 2*Math.abs(kappa())*cross(xIn, yIn)/xi()
+ );
+ }
+
+ private double aTwo(double xIn, double yIn){
+ return ((double)
+ 2*kappa()*( alpha()*xOn(xIn, yIn) + beta()*yOn(xIn, yIn) )/xi() + xi()
+ );
+ }
+
+ private double sPerp(double xIn, double yIn){
+
+ double daOne = aOne(xIn, yIn);
+ double daTwo = aTwo(xIn, yIn);
+
+ if ( daOne >= 0 && daOne < 0.3 ){
+ return ((double)
+ cross(xIn, yIn)
+ *(1+daOne*daOne*(1/6 + daOne*daOne*(3/40)))
+ /xi()
+ );
+ } else {
+ double temp = Math.atan2(daOne,daTwo);
+ if ( temp < 0 ) temp += 2 * Math.PI;
+ return ((double)
+ temp / ( 2*Math.abs(kappa()) )
+ );
+ }
+ }
+
+ private double lambda(){
+ return ((double)
+ (Z[1]-Z[2])/(sPerp(X[1], Y[1])-sPerp(X[2], Y[2]))
+ );
+ }
+ private double z0(){
+ return ((double)
+ Z[1] - lambda()*sPerp(X[1], Y[1])
+ );
+ }
+
+
+//
+// Field Summary
+//
+ private int numberOfHits;
+ private double[][] Hits = new double[100][3];
+ private double[][] errorHits = new double[100][3];
+
+ private MIPCluster thisMIPObject;
+ private CalorimeterCell thisCell;
+
+ private TrackParameters tP = null;
+
+ private double[] X;
+ private double[] Y;
+ private double[] Z;
+
+}
lcd/hep/lcd/recon/cluster/MIPMatching
diff -N ThreeDChi2.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ThreeDChi2.java 10 Aug 2005 17:38:18 -0000 1.1
@@ -0,0 +1,44 @@
+package hep.lcd.recon.cluster.MIPMatching;
+
+import java.io.*;
+import java.util.*;
+
+/**
+ 3D Chi2 for MIP Matching
+
+ @version 1.0
+ @author [log in to unmask]
+
+ Modification Log:
+ -- 08/08/2004 Version 1.0
+*/
+
+public class ThreeDChi2 extends AbstractChi2{
+
+ public double getChi2(){
+ //
+ // d0 contribution
+ //
+ double d0chi2 = (tp1.d0()-tp2.d0())*(tp1.d0()-tp2.d0()) / (ed0*ed0);
+ //
+ // phi0 contribution
+ //
+ double phi0chi2 = (tp1.phi0()-tp2.phi0())*(tp1.phi0()-tp2.phi0()) / (ephi0*ephi0);
+ //
+ // kappa contribution
+ //
+ double kappachi2 = (tp1.kappa()-tp2.kappa())*(tp1.kappa()-tp2.kappa()) / (ekappa*ekappa);
+ //
+ // tanl contribution
+ //
+ double tanlchi2 = (tp1.tanl()-tp2.tanl())*(tp1.tanl()-tp2.tanl()) / (etanl*etanl);
+ //
+ // z0 contribution
+ //
+ double z0chi2 = (tp1.z0()-tp2.z0())*(tp1.z0()-tp2.z0()) / (ez0*ez0);
+
+ double chi2 = d0chi2 + phi0chi2 + kappachi2 + tanlchi2 + z0chi2;
+ return ((double) chi2);
+ };
+
+}
lcd/hep/lcd/recon/cluster/MIPMatching
diff -N TrackParameters.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackParameters.java 10 Aug 2005 17:38:18 -0000 1.1
@@ -0,0 +1,269 @@
+package hep.lcd.recon.cluster.MIPMatching;
+/**
+ Track Parameters in (kappa, d0, phi0, z0, tanl)-representation.
+
+ @author Wolfgang F. Mader ([log in to unmask])
+ @version 1.0
+
+ Modification Log:
+ -- 07/22/2004 Version 1.0
+*/
+
+//public class TPkdpzt extends TrackParameters{
+public class TrackParameters{
+//
+// Constructor Summary
+//
+ public TrackParameters(){
+ setKappa(0.001);
+ setD0(0.0001);
+ setPhi0(0.0001);
+ setZ0(0.0001);
+ setTanl(0.0001);
+ }
+
+ public TrackParameters(double k, double d, double p, double z, double tl){
+ setKappa(k);
+ setD0(d);
+ setPhi0(p);
+ setZ0(z);
+ setTanl(tl);
+ }
+
+ public TrackParameters(double[] d){
+ this(d[0], d[1], d[2], d[3], d[4]);
+ }
+//
+// Method Summary
+//
+//..Provide OPAL representation
+ private void setKappa(double d){
+ trackParameter[0] = d;
+ }
+ private void setD0(double d){
+ trackParameter[1] = d;
+ }
+ private void setPhi0(double d){
+ trackParameter[2] = d;
+ }
+ private void setZ0(double d){
+ trackParameter[3] = d;
+ }
+ private void setTanl(double d){
+ trackParameter[4] = d;
+ }
+
+ public double kappa(){
+ return trackParameter[0];
+ }
+ public double d0(){
+ return trackParameter[1];
+ }
+ public double phi0(){
+ return trackParameter[2];
+ }
+ public double z0(){
+ return trackParameter[3];
+ }
+ public double tanl(){
+ return trackParameter[4];
+ }
+
+//.. Provide alpha beta gamma representation
+ public double alpha(){
+ double phi0 = trackParameter[2];
+ return((double)
+ Xid()*Math.sin(phi0)
+ );
+ }
+ public double beta(){
+ double phi0 = trackParameter[2];
+ return ((double)
+ Xid()*Math.cos(phi0)
+ );
+ }
+ public double gamma(){
+ double d0 = trackParameter[1];
+ return ((double)
+ d0 * (1.0 + akad())
+ );
+ }
+ private double akad(){
+ double kappa = trackParameter[0];
+ double d0 = trackParameter[1];
+ return ((double)
+ kappa * d0
+ );
+ }
+ private double Xid(){
+ return((double)
+ (1.0 + 2.0*akad())
+ );
+ }
+
+//.. Provide W-Presentation
+ public double pX(){
+ double kappa = Math.abs(kappa());
+ double phi0 = phi0();
+
+ return ((double)
+ (A()/kappa)*Math.cos(phi0)
+ );
+ }
+ public double pY(){
+ double kappa = Math.abs(kappa());
+ double phi0 = phi0();
+ return ((double)
+ (A()/kappa)*Math.sin(phi0)
+ );
+ }
+ public double pZ(){
+ double kappa = Math.abs(kappa());
+ double tanl = tanl();
+
+ return ((double)
+ (A()/kappa)*tanl
+ );
+ }
+ public double X(){
+ double d0 = d0();
+ double phi0 = phi0();
+ return ((double)
+ -d0 * Math.sin(phi0)
+ );
+ }
+ public double Y(){
+ double d0 = d0();
+ double phi0 = phi0();
+ return ((double)
+ d0 * Math.cos(phi0)
+ );
+ }
+ public double Z(){
+ double z0 = z0();
+ return ((double) z0);
+ }
+
+ private double A(){
+ return ((double)
+ 0.00015 * 50.
+ );
+ }
+
+ public void moveTrackParameters(double[] x){
+ moveTrackParameters(x[0], x[1], x[2]);
+ }
+
+ public void moveTrackParameters(double x, double y, double z){
+ double[] temp = getTrackParametersAt(x, y, z);
+
+ setKappa(temp[0]);
+ setD0(temp[1]);
+ setPhi0(temp[2]);
+ setZ0(temp[3]);
+ setTanl(temp[4]);
+ }
+
+ public double[] getParameters(){
+ return trackParameter;
+ }
+
+ public double[] getTrackParametersAt(double x, double y, double z){
+
+ double[] parameters = new double[] {-999.,-999.,-999.,-999.,-999};
+
+ double AKAIN = kappa();
+ double PHI0IN = phi0();
+ double D0IN = d0();
+ double TANLIN = tanl();
+ double Z0IN = z0();
+
+ double ALFAI = alpha();
+ double BETAI = beta();
+ double GAMMAI = gamma();
+ double XSII = Math.sqrt(1.0+4*AKAIN*GAMMAI);
+
+ parameters[0] = AKAIN;
+
+ double GAMMAO = AKAIN*(Math.pow(x,2)+Math.pow(y,2))+ALFAI*x+BETAI*y+GAMMAI;
+ double DELH = GAMMAO*(1. - AKAIN*GAMMAO);
+//
+// use an exact calculation if we need to...
+//
+ double AKADEL = 4.*AKAIN*GAMMAO;
+ if (Math.abs(DELH*AKAIN) > 0.01 && AKADEL > -1.){
+ DELH = 0.5/AKAIN*(Math.sqrt(1.+ AKADEL) - 1.);
+ }
+
+ double XXON = (x - DELH*ALFAI)/(1.+2.*AKAIN*DELH);
+ double YYON = (y - DELH*BETAI)/(1.+2.*AKAIN*DELH);
+ double XON = XXON;
+ double YON = YYON;
+
+ double YYTMP = ( 2.*AKAIN)*XXON + ALFAI;
+ double XXTMP = (-2.*AKAIN)*YYON - BETAI;
+ double PHID = Math.atan2(YYTMP,XXTMP);
+
+ parameters[1] = DELH;
+
+ if ( PHID < 0.0 ) PHID += 2.*Math.PI;
+
+ parameters[2] = PHID;
+
+ double XSIO = 1. + DELH*2.*AKAIN;
+//
+// Calculate arc length (in x-y plane)
+// (choose output range from -2*pi*R to +2*pi*R)
+//
+ double DELFI = (PHID - PHI0IN) % (2.*Math.PI);
+ if ( Math.abs(DELFI) > Math.PI && DELFI >= 0){
+ DELFI = DELFI - 2.*Math.PI;
+ } else if ( Math.abs(DELFI) > Math.PI && DELFI < 0){
+ DELFI = DELFI + 2.*Math.PI;
+ }
+
+
+ double SS = -999.;
+ if ( Math.abs(AKAIN) > 0.000001 ){
+ SS = DELFI/(2.*AKAIN);
+ } else {
+// straight line limit (akain -> 0 ; xsii -> 1 ; delfi -> 0 )
+ SS = Math.cos(PHID)*XXON + Math.sin(PHID)*YYON;
+ }
+
+ if ( Math.abs(SS) < 0.00000001) SS = 0.;
+//
+// Now update Z parameters
+//
+ double ZON = Z0IN + SS*TANLIN;
+ parameters[3] = ZON - z;
+ parameters[4] = TANLIN;
+
+ return parameters;
+ }
+
+ public void updateParameters(double k, double d, double p, double z, double tl){
+ setKappa(k);
+ setD0(d);
+ setPhi0(p);
+ setZ0(z);
+ setTanl(tl);
+ }
+ public void updateParameters(double[] d){
+ setKappa(d[0]);
+ setD0(d[1]);
+ setPhi0(d[2]);
+ setZ0(d[3]);
+ setTanl(d[4]);
+ }
+
+ public double[][] getErrorMatrix(){
+ return errorMatrix;
+ }
+//
+// Field Summary
+//
+ double[] referencePoint = new double[3];
+ double[] trackParameter = new double[5];
+ double[][] errorMatrix = new double[5][5];
+}
CVSspam 0.2.8