lcsim/src/org/lcsim/contrib/RobKutschke/TRFSelfTest/hits
diff -N RKHit.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ RKHit.java 8 Jun 2009 06:04:24 -0000 1.1
@@ -0,0 +1,267 @@
+package org.lcsim.contrib.RobKutschke.TRFSelfTest.hits;
+
+import java.util.List;
+import java.util.Random;
+
+
+import org.lcsim.recon.tracking.trfbase.PropStat;
+import org.lcsim.recon.tracking.trfbase.VTrack;
+import org.lcsim.recon.tracking.trfbase.Hit;
+import org.lcsim.recon.tracking.trfbase.TrackError;
+import org.lcsim.recon.tracking.trfbase.ETrack;
+
+import org.lcsim.recon.tracking.trfcyl.ClusCylPhi;
+import org.lcsim.recon.tracking.trfcyl.ClusCylPhiZ2D;
+import org.lcsim.recon.tracking.trfcyl.HitCylPhi;
+import org.lcsim.recon.tracking.trfcyl.HitCylPhiZ2D;
+import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
+
+import org.lcsim.recon.tracking.trfzp.ClusZPlane1;
+import org.lcsim.recon.tracking.trfzp.ClusZPlane2;
+import org.lcsim.recon.tracking.trfzp.HitZPlane1;
+import org.lcsim.recon.tracking.trfzp.HitZPlane2;
+import org.lcsim.recon.tracking.trfzp.SurfZPlane;
+
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.generator.RKTrack;
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.geom.RKGeom;
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.geom.RKSurf;
+
+import org.lcsim.recon.tracking.trfutil.TRFMath;
+
+import org.lcsim.util.aida.AIDA;
+/**
+ *
+ * Interface between RKGeom and the TRF hit classes.
+ * The main purpose of this class is to produce the
+ * various sorts of TRF hits.
+ *
+ *@author $Author: kutschke $
+ *@version $Id: RKHit.java,v 1.1 2009/06/08 06:04:24 kutschke Exp $
+ *
+ * Date $Date: 2009/06/08 06:04:24 $
+ *
+ */
+
+public class RKHit{
+
+ static private Random ran = new Random();
+ static private Random ranz = new Random();
+
+ private AIDA aida = AIDA.defaultInstance();
+
+
+ // Resolutions in cm ( ie TRF units ).
+ //static private double respixel = 0.0025;
+ //static private double resstrip = 0.0050;
+ //static private double zres_strip = 2.89; // 10 cm /sqrt(12).
+ //static private double zres_strip = 1.00; // see if a smaller number helps.
+
+ private RKTrack track;
+ private RKSurf surf;
+ private VTrack vt;
+ private PropStat ps;
+ private double path;
+ private double psi;
+ private Hit ghit = null;
+
+ // Dimension of the readout
+ private int rodim = 0;
+
+ // Number of stereo measurements in this device.
+ private int zdim = 0;
+
+ public RKHit( RKTrack track, RKSurf surf, VTrack vt, PropStat ps, double path){
+ this.track = track;
+ this.surf = surf;
+ this.vt = new VTrack(vt);
+ this.ps = new PropStat(ps);
+ this.path = path;
+ psi = path*track.sz()/track.rho();
+
+ rodim = surf.rodim;
+ if ( surf.getType() == RKSurf.type_zdisc ){
+ zdim = 1;
+ } else if ( surf.getType() == RKSurf.type_tube ){
+ if ( rodim == 2 ) {
+ zdim=1;
+ }
+ }
+
+ }
+
+
+ public RKTrack getRKTrack() { return track; }
+ public RKSurf getSurface() { return surf; }
+ public VTrack getVTrack() { return vt; }
+ public PropStat getPropStat(){ return ps; }
+ public double getPath() { return path; }
+ public double getPsi() { return psi; }
+ public int getRodim() { return rodim; }
+ public int getZdim() { return zdim; }
+
+ static public void setSeed( long seed ){
+ ran.setSeed(seed);
+ ranz.setSeed(seed+1234578);
+ }
+
+ public Hit MakeHit(){
+ if ( ghit != null ) return ghit;
+
+ if ( surf.getType() == RKSurf.type_tube ){
+ if ( surf.rodim == 0 ){
+ return null;
+ } else if ( surf.rodim == 1 ){
+ return MakeCylPhi( surf.getResolution0() );
+ } else if ( surf.rodim == 2 ){
+ return MakeCylPhiZ2D( surf.getResolution0(), surf.getResolution1() );
+ }
+ } else if ( surf.getType() == RKSurf.type_zdisc ){
+ if ( surf.rodim == 0 ){
+ return null;
+ } else if ( surf.rodim == 1 ){
+ return MakeZPlane1D( surf.getResolution0() );
+ } else if ( surf.rodim == 2 ){
+ return MakeZPlane2D( surf.getResolution0(), surf.getResolution1() );
+ }
+ }
+ return null;
+ }
+
+ private Hit MakeCylPhi( double res){
+ SurfCylinder s = (SurfCylinder)vt.surface();
+ double r = s.radius();
+
+ double sig0 = res/r;
+
+ double m0 = vt.vector(0) + sig0*ran.nextGaussian();
+ //double m0 = vt.vector(0);
+
+ double stereo = 0.;
+ int mcid = 0;
+ ClusCylPhi cluster = new ClusCylPhi( r, m0, sig0, mcid );
+ TrackError verr = new TrackError();
+ verr.set(0,0,10.);
+ verr.set(1,1,100.);
+ verr.set(2,2,2.);
+ verr.set(3,3,10.);
+ verr.set(4,4,10.);
+
+ ETrack et = new ETrack(s.newPureSurface(), vt.vector(), verr );
+ List ghits = cluster.predict(et,cluster);
+ ghit = (HitCylPhi)ghits.get(0);
+ return ghit;
+ }
+
+ private Hit MakeCylPhiZ2D( double sigpixel, double sigz ){
+
+ SurfCylinder s = (SurfCylinder)vt.surface();
+ double r = s.radius();
+
+ double sig0 = sigpixel/r;
+ double sig1 = sigz;
+
+ double m0 = vt.vector(0) + sig0*ran.nextGaussian();
+
+ //double m1 = vt.vector(1) + sig1*ran.nextGaussian();
+
+ // Hack to decouple r-phi and z fits.
+ double m1;
+ if ( r > 100 ){
+ m1 = vt.vector(1) + sig1*ranz.nextGaussian();
+ }else{
+ m1 = vt.vector(1) + sig1*ran.nextGaussian();
+ }
+
+ //double m0 = vt.vector(0);
+ //double m1 = vt.vector(1);
+
+ double stereo = 0.;
+ int mcid = 0;
+ ClusCylPhiZ2D cluster = new ClusCylPhiZ2D( r, m0, sig0, m1, sig1, stereo, mcid );
+ TrackError verr = new TrackError();
+ verr.set(0,0,10.);
+ verr.set(1,1,100.);
+ verr.set(2,2,2.);
+ verr.set(3,3,10.);
+ verr.set(4,4,10.);
+
+ ETrack et = new ETrack(s.newPureSurface(), vt.vector(), verr );
+ List ghits = cluster.predict(et,cluster);
+ ghit = (HitCylPhiZ2D)ghits.get(0);
+ return ghit;
+ }
+
+ private Hit MakeZPlane1D( double res){
+ SurfZPlane s = (SurfZPlane)vt.surface();
+ double z = s.z();
+
+ // Measurement error.
+ double sig0 = res;
+
+ // Measurement direction.
+ double wx = 1.;
+ double wy = 0.;
+ if ( surf.ixy == RKSurf.ixy_y ){
+ wx = 0.;
+ wy = 1.;
+ }
+ double m0 = vt.vector(0)*wx + vt.vector(1)*wy;
+ m0 += sig0*ran.nextGaussian();
+
+ int mcid = 0;
+ ClusZPlane1 cluster = new ClusZPlane1( z, wx, wy, m0, sig0, mcid );
+ TrackError verr = new TrackError();
+ verr.set(0,0,10.);
+ verr.set(1,1,10.);
+ verr.set(2,2,10.);
+ verr.set(3,3,10.);
+ verr.set(4,4,10.);
+
+ ETrack et = new ETrack(s.newPureSurface(), vt.vector(), verr );
+ List ghits = cluster.predict(et,cluster);
+ ghit = (HitZPlane1)ghits.get(0);
+ return ghit;
+
+ }
+
+ private Hit MakeZPlane2D( double res0, double res1){
+ SurfZPlane s = (SurfZPlane)vt.surface();
+ double z = s.z();
+
+ double sig0 = res0;
+ double sig1 = res1;
+
+ double m0 = vt.vector(0) + sig0*ran.nextGaussian();
+ double m1 = vt.vector(1) + sig1*ran.nextGaussian();
+ //double m0 = vt.vector(0);
+ //double m1 = vt.vector(1);
+
+ double dxdy = 0.;
+ int mcid = 0;
+ ClusZPlane2 cluster = new ClusZPlane2( z, m0, m1, sig0*sig0, sig1*sig1, dxdy, mcid );
+ TrackError verr = new TrackError();
+ verr.set(0,0,10.);
+ verr.set(1,1,10.);
+ verr.set(2,2,10.);
+ verr.set(3,3,10.);
+ verr.set(4,4,10.);
+
+ ETrack et = new ETrack(s.newPureSurface(), vt.vector(), verr );
+ List ghits = cluster.predict(et,cluster);
+ ghit = (HitZPlane2)ghits.get(0);
+ return ghit;
+
+ }
+
+ private Hit MakeNullCyl( ){
+ SurfCylinder s = (SurfCylinder)vt.surface();
+ double r = s.radius();
+
+ ghit = new HitNull();
+
+ return ghit;
+ }
+
+
+
+}
lcsim/src/org/lcsim/contrib/RobKutschke/TRFSelfTest/hits
diff -N RKHitList.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ RKHitList.java 8 Jun 2009 06:04:24 -0000 1.1
@@ -0,0 +1,288 @@
+package org.lcsim.contrib.RobKutschke.TRFSelfTest.hits;
+
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.geom.RKSurf;
+
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Collections;
+import java.util.Comparator;
+
+/**
+ *
+ * Holds a list of RKHits. Can return the hits sorted in different
+ * orders. Can return properties of the ensemble of hits.
+ *
+ *@author $Author: kutschke $
+ *@version $Id: RKHitList.java,v 1.1 2009/06/08 06:04:24 kutschke Exp $
+ *
+ * Date $Date: 2009/06/08 06:04:24 $
+ *
+ */
+
+public class RKHitList {
+
+ // Original list.
+ private List<RKHit> list = new ArrayList<RKHit>();
+
+ // Sorted copies of the lists.
+ private List<RKHit> forward = null;
+ private List<RKHit> backward = null;
+
+ private int nmeas = 0;
+ private int nzmeas = 0;
+
+ private int ncyl = 0;
+ private int nzp = 0;
+
+ RKHitList (){
+ }
+
+ public void addHit( RKHit h ){
+ list.add(h);
+ nmeas += h.getRodim();
+ nzmeas += h.getZdim();
+
+ int type = h.getSurface().getType();
+ if ( type == RKSurf.type_tube ){
+ ncyl += 1;
+ } else if ( type == RKSurf.type_zdisc ){
+ nzp += 1;
+ }
+
+ // Any existing sorts are no longer valid.
+ forward = null;
+ backward = null;
+ }
+
+ public RKHit get( int i){
+ return list.get(i);
+ }
+
+ public List<RKHit> getForward(){
+ if( forward != null ) return forward;
+ forward = new ArrayList<RKHit>(list);
+ Collections.sort( forward, new CompareRKHit<RKHit>(+1) );
+ return forward;
+ }
+
+ public List<RKHit> getBackward(){
+ if ( backward != null ) return backward;
+ backward = new ArrayList<RKHit>(list);
+ Collections.sort( backward, new CompareRKHit<RKHit>(-1) );
+ return backward;
+ }
+
+ public int nHits(){ return list.size(); }
+ public int nDof(){ return nDof(5); }
+ public int nDof( int n){ return nmeas-n; }
+
+ public int nZ(){ return nzmeas;}
+
+ public int nCyl(){ return ncyl;}
+ public int nZp() { return nzp;}
+
+ // Does the track have enough measurements to be fittable.
+ // This is an arbitrary definition.
+ public boolean fitable(){
+ if ( nDof() < 1 ) return false;
+ if ( nzmeas < 3 ) return false;
+ return true;
+ }
+
+
+ public RKHit outerMostHit(){
+ for ( RKHit h : getBackward() ){
+ if ( h.getRodim() > 0 ) return h;
+ }
+ return null;
+ }
+
+ public RKHit innerMostHit(){
+ for ( RKHit h : getForward() ){
+ if ( h.getRodim() > 0 ) return h;
+ }
+ return null;
+ }
+
+ public int nLeadingStrips(){
+ int n = 0;
+ int nLeadingStrips = 0;
+ int firstZType = 0;
+ boolean firstZ = false;
+ for ( RKHit h: getBackward() ){
+ RKSurf s = h.getSurface();
+ int rodim = h.getRodim();
+ int zdim = h.getZdim();
+ if ( zdim > 0 ){
+ firstZ = true;
+ if ( s.getType() == RKSurf.type_tube ){
+ firstZType = 1;
+ } else{
+ firstZType = 2;
+ }
+ break;
+ }
+ if ( rodim == 1 &
+ s.getType() == RKSurf.type_tube ){
+ ++nLeadingStrips;
+ }
+
+ }
+ return nLeadingStrips;
+ }
+
+
+ public int firstZType(){
+ int type = 0;
+ for ( RKHit h: getBackward() ){
+ RKSurf s = h.getSurface();
+ int zdim = h.getZdim();
+ if ( zdim > 0 ){
+ if ( s.getType() == RKSurf.type_tube ){
+ type = 1;
+ } else{
+ type = 2;
+ }
+ break;
+ }
+ }
+ return type;
+ }
+
+ public void PrintBackward(){
+ for ( RKHit h: getBackward() ){
+ RKSurf s = h.getSurface();
+ s.Print();
+ /*
+ int zdim = h.getZdim();
+ int rodim = h.getRodim();
+ String ctype = "";
+ double zc = 0.;
+ if ( rodim > 0 ){
+ if ( s.getType() == RKSurf.type_tube ){
+ ctype = "Barrel";
+ } else{
+ ctype = "Endcap";
+ zc = s.zc;
+ }
+ System.out.printf ( "Hit: %6s %4d %4d %9.2f %9.2f\n",
+ ctype,
+ rodim,
+ s.ixy,
+ h.getPath(),
+ zc
+ );
+
+ } else{
+ s.Print();
+ }
+ */
+ }
+ }
+
+ public double deltaZ(){
+ RKHit h0 = getBackward().get(0);
+ RKSurf s0 = h0.getSurface();
+ int ro0 = h0.getRodim();
+
+ // Only consider tracks that start with zdisc strips.
+ if ( s0.getType() != RKSurf.type_zdisc ) return 0.;
+ if ( ro0 != 1 ) return 0.;
+
+ // z of first z measuring surface.
+ double z0 = s0.zc;
+
+ // z of first z measuring surface with enough lever arm to be useful.
+ double z1 = z0;
+
+ int n = 0;
+ for ( RKHit h: getBackward() ){
+
+ // Start looking at the third measurement.
+ // ie skip second plane if it is part of a doublet.
+ if ( ++n < 3 ) continue;
+
+ RKSurf s = h.getSurface();
+ if ( s.getType() == RKSurf.type_zdisc ){
+ // Accept any z disc.
+ z1 = s.zc;
+ break;
+ } else if ( h.getRodim() > 1 ) {
+ // Accept barrel pixels.
+ z1 = h.getVTrack().vector(1);
+ break;
+ }
+
+ }
+
+ return Math.abs(z0-z1);
+ }
+
+ // Look for patterns of z hits that can cause problems.
+ public int Pattern(){
+
+ // Number of endcap hits before the first barrel hit.
+ int nec = 0;
+ for ( RKHit h: getBackward() ){
+ RKSurf s = h.getSurface();
+ if ( s.getType() == RKSurf.type_zdisc ){
+ ++nec;
+ }else{
+ break;
+ }
+ }
+
+ // The first hit is a barrel hit.
+ if ( nec == 0 ) return 0;
+
+ // There is are hits from at least 2 z stations so at least
+ // one slope will be well defined.
+ if ( nec > 2 ) return 0;
+
+ int nbar = 0;
+ int nskip = 0;
+ for ( RKHit h: getBackward() ){
+ if ( ++nskip <= nec ) continue;
+
+ RKSurf s = h.getSurface();
+ if ( s.getType() == RKSurf.type_tube ){
+ ++nbar;
+ }else{
+ break;
+ }
+ }
+
+ // Number of barrel hits between the first z hit and
+ return nbar;
+
+ }
+
+ // ??? This should be in the RKHits definition, not in this file????
+ //
+ // Class to compare two RKHits, sorted by path length.
+ // If the argument is positive, then the sort is done in increasing order.
+ // If the argument is negative, then the sort is done in decreasing order.
+ private class CompareRKHit<T> implements Comparator<T>{
+ private int sign = 1;
+ public CompareRKHit(){
+ }
+ public CompareRKHit(int sign){
+ this.sign = sign;
+ }
+ public int compare ( Object o1, Object o2 ){
+ RKHit h1 = (RKHit)o1;
+ RKHit h2 = (RKHit)o2;
+ double l1 = h1.getPath();
+ double l2 = h2.getPath();
+ if ( l1 == l2 ) {
+ return 0;
+ } else if ( l1 > l2 ) {
+ return +1*sign;
+ } else{
+ return -1*sign;
+ }
+ }
+ }
+
+
+}
lcsim/src/org/lcsim/contrib/RobKutschke/TRFSelfTest/hits
diff -N RKMakeHits.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ RKMakeHits.java 8 Jun 2009 06:04:24 -0000 1.1
@@ -0,0 +1,645 @@
+package org.lcsim.contrib.RobKutschke.TRFSelfTest.hits;
+
+import java.util.List;
+import java.util.Random;
+import Jama.Matrix;
+
+import org.lcsim.util.aida.AIDA;
+
+import org.lcsim.recon.tracking.trfbase.TrackError;
+import org.lcsim.recon.tracking.trfbase.ETrack;
+import org.lcsim.recon.tracking.trfbase.Hit;
+import org.lcsim.recon.tracking.trfbase.VTrack;
+import org.lcsim.recon.tracking.trfbase.PropStat;
+import org.lcsim.recon.tracking.trfbase.PropDir;
+import org.lcsim.recon.tracking.trfbase.Propagator;
+
+import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
+import org.lcsim.recon.tracking.trfzp.SurfZPlane;
+
+import org.lcsim.recon.tracking.trfcyl.ThinCylMs;
+import org.lcsim.recon.tracking.trfcyl.ThinCylMsSim;
+//import org.lcsim.recon.tracking.trfcyl.CylEloss;
+//import org.lcsim.recon.tracking.trfcyl.CylElossSim;
+
+
+import org.lcsim.contrib.RobKutschke.TRF.trfcyl.CylEloss;
+import org.lcsim.contrib.RobKutschke.TRF.trfcyl.CylElossSim;
+import org.lcsim.contrib.RobKutschke.TRF.trfzp.ThinZPlaneMs;
+import org.lcsim.contrib.RobKutschke.TRF.trfzp.ThinZPlaneMsSim;
+import org.lcsim.contrib.RobKutschke.TRF.trfeloss.RKDeDxFixed;
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.geom.RKGeom;
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.geom.RKSurf;
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.util.RKDebug;
+
+//import org.lcsim.recon.tracking.trfeloss.DeDxBethe;
+import org.lcsim.recon.tracking.trfeloss.DeDx;
+import org.lcsim.recon.tracking.trfeloss.DeDxFixed;
+
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.generator.toTRF;
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.generator.RKTrack;
+
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.util.VTUtil;
+
+import org.lcsim.contrib.RobKutschke.ToyConfig.*;
+
+/**
+ * Propagate a track through a detector and generate hits.
+ * Include multiple scattering and eloss, or not, as controlled by
+ * the configuration file.
+ *
+ *@author $Author: kutschke $
+ *@version $Id: RKMakeHits.java,v 1.1 2009/06/08 06:04:24 kutschke Exp $
+ *
+ * Date $Date: 2009/06/08 06:04:24 $
+ *
+ */
+
+
+public class RKMakeHits {
+
+ // Copies of arguments to c'tor.
+ // The geometry summary extracted from the GeomConverter classes.
+ RKGeom rkgeom = null;
+
+ // A random number generator to use.
+ Random random;
+
+ private AIDA aida = AIDA.defaultInstance();
+
+ // TRF wants distances in cm, not mm.
+ private double mmTocm = 0.1;
+
+ // Thicknesses in radiation lengths of the 3 classes of cylindrical surfaces.
+ private double scatThick0 = 0.006136;
+ private double scatThick1 = 0.000916;
+ private double scatThick2 = 0.013747;
+
+ // Simulators for scattering at the 3 classes of cylindrical surfaces
+ // and 2 classes of z surfaces.
+ private ThinCylMsSim scatSim0 = null;
+ private ThinCylMsSim scatSim1 = null;
+ private ThinCylMsSim scatSim2 = null;
+ private ThinZPlaneMsSim zscatSim1 = null;
+ private ThinZPlaneMsSim zscatSim2 = null;
+
+ // Simulators for energy loss at the 3 classes of cylindrical surfaces
+ // and 2 classes of z surfaces.
+ private CylElossSim celoss0 = null;
+ private CylElossSim celoss1 = null;
+ private CylElossSim celoss2 = null;
+ private CylElossSim zeloss1 = null;
+ private CylElossSim zeloss2 = null;
+
+ // The run time configuration system.
+ ToyConfig config;
+
+ // Quantities taken from the run time configuration system.
+ private boolean doms = false;
+ private boolean doeloss = false;
+ private double msfac = 1.0;
+ private double dedxscale = 7.5;
+ private double dedxsigma = 0.00;
+ private boolean doHitCheck = false;
+
+ public RKMakeHits( RKGeom rkgeom, Random random ){
+ this.rkgeom = rkgeom;
+ this.random = random;
+
+ // Copy information from the run time configuration system.
+ try{
+ ToyConfig config = ToyConfig.getInstance();
+ doms = config.getBoolean("DoMS");
+ doeloss = config.getBoolean("DoELoss");
+ msfac = config.getDouble("MsFac");
+ dedxscale = config.getDouble("dEdXScale");
+ dedxsigma = config.getDouble("dEdXSigma");
+ doHitCheck = config.getBoolean("doHitCheck");
+
+ } catch (ToyConfigException e){
+ System.out.println (e.getMessage() );
+ System.out.println ("Stopping now." );
+ System.exit(-1);
+ }
+ System.out.println ("RKMakeHits dedxscale: " + dedxscale );
+
+ // Seed for track random number generator.
+ long seed = -398783512;
+
+ // An increment used to generate independent streams of random
+ // numbers - I have no proof that this really works.
+ long seedIncrement = 876633217;
+
+ CylEloss.SetNewModel(true);
+
+ scatSim0 = new ThinCylMsSim(scatThick0*msfac);
+ scatSim1 = new ThinCylMsSim(scatThick1*msfac);
+ scatSim2 = new ThinCylMsSim(scatThick2*msfac);
+
+ zscatSim1 = new ThinZPlaneMsSim(scatThick1*msfac);
+ zscatSim2 = new ThinZPlaneMsSim(scatThick2*msfac);
+
+ // Seed the random number generators in the Siminteractors.
+ seed += seedIncrement;
+ scatSim0.setSeed(seed);
+ seed += seedIncrement;
+ scatSim1.setSeed(seed);
+ seed += seedIncrement;
+ scatSim2.setSeed(seed);
+ seed += seedIncrement;
+ zscatSim1.setSeed(seed);
+ seed += seedIncrement;
+ zscatSim2.setSeed(seed);
+
+ double densityBe = rkgeom.getCylinders().get(0).density;
+ double densitySi = rkgeom.getCylinders().get(1).density;
+ RKDeDxFixed dedxBe = new RKDeDxFixed(densityBe,dedxscale, dedxsigma);
+ RKDeDxFixed dedxSi = new RKDeDxFixed(densitySi,dedxscale, dedxsigma);
+
+ double thickbeampipe = rkgeom.getCylinders().get(0).thick;
+ double thickpixel = rkgeom.getCylinders().get(1).thick;
+ double thickstrip = rkgeom.getCylinders().get(6).thick;
+
+ celoss0 = new CylElossSim( thickbeampipe, dedxBe, random );
+ celoss1 = new CylElossSim( thickpixel, dedxSi, random );
+ celoss2 = new CylElossSim( thickstrip, dedxSi, random );
+
+ zeloss1 = new CylElossSim( thickpixel, dedxSi, random );
+ zeloss2 = new CylElossSim( thickstrip, dedxSi, random );
+ System.out.println( "ELoss thicknesses: " + thickbeampipe + " " + thickpixel + " " + thickstrip );
+ System.out.println( "Densities: " + densityBe + " " + densitySi );
+
+
+ // This is probably not necessary - check later.
+ seed += seedIncrement;
+ celoss0.setSeed(seed);
+ seed += seedIncrement;
+ celoss1.setSeed(seed);
+ seed += seedIncrement;
+ celoss2.setSeed(seed);
+ seed += seedIncrement;
+ zeloss1.setSeed(seed);
+ seed += seedIncrement;
+ zeloss2.setSeed(seed);
+
+ }
+
+ // Generate perfect hits on the outward going arc only.
+ // Do not apply energy loss or multiple scattering.
+ public RKHitList GenerateOneArcHits( RKTrack rktrk ){
+
+ // Convert to TRF style track.
+ toTRF trftrk = new toTRF(rktrk);
+ VTrack vtdca = trftrk.atDcaZaxis();
+ VTrack vtz = trftrk.atZPlane(rktrk.z0());
+
+ // This code only looks at the first arc of a track.
+ double pathlimit = trftrk.halfarc();
+ int nmeas = 0;
+ int nzmeas = 0;
+ int ncmeas = 0;
+ RKHitList rkl = new RKHitList();
+ Propagator prop = rkgeom.newPropagator();
+ double sum = 0.;
+ for ( RKSurf s: rkgeom.getCylinders() ){
+ PropStat ps = prop.vecDirProp( vtdca, s.getCylinder(), PropDir.FORWARD);
+ if ( ps.success() ){
+ sum += ps.pathDistance();
+ if ( sum > pathlimit ) break;
+ if ( s.inBounds(vtdca.vector(1)) ){
+ rkl.addHit( new RKHit( rktrk, s, vtdca, ps, sum ) );
+ aida.cloud2D( "HitGenNoScat/C r vs z",-1).fill( vtdca.vector(1), s.radius*mmTocm );
+ aida.cloud2D( "HitGenNoScat/Both r vs z",-1).fill( vtdca.vector(1), s.radius*mmTocm );
+ double r = s.radius*mmTocm;
+ /*
+ System.out.printf ("Adding Cyl: %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n",
+ r*Math.cos(vtdca.vector(0)),
+ r*Math.sin(vtdca.vector(0)),
+ r,
+ vtdca.vector(1),
+ sum,
+ vtdca.vector(3)
+ );
+ */
+ nmeas += s.rodim;
+ ncmeas += s.rodim;
+ }
+ }
+ }
+ aida.cloud2D("HitGenNoScat/Path length vs cz",-1).fill(rktrk.cz(),sum);
+ aida.cloud1D("HitGenNoScat/Sum").fill(sum);
+ aida.cloud2D("HitGenNoScat/Hits vs cz",-1).fill( rktrk.cz(), ncmeas );
+
+ sum = 0;
+ for ( RKSurf s: rkgeom.getZ( rktrk.z0(), rktrk.cz() ) ){
+ PropStat ps = prop.vecDirProp( vtz, s.getZPlane(), PropDir.FORWARD );
+ if ( ps.success() ){
+ sum += ps.pathDistance();
+ if ( sum > pathlimit ) break;
+ double r = Math.sqrt(vtz.vector(0)*vtz.vector(0) + vtz.vector(1)*vtz.vector(1) );
+ if ( s.inBounds(r) ){
+ rkl.addHit( new RKHit( rktrk, s, vtz, ps, sum ) );
+ aida.cloud2D( "HitGenNoScat/Z r vs z",-1).fill( s.zc*mmTocm, r );
+ aida.cloud2D( "HitGenNoScat/Both r vs z",-1).fill( s.zc*mmTocm, r );
+ nmeas += s.rodim;
+ nzmeas += s.rodim;
+ /*
+ System.out.printf ("Adding Z: %10.4f %10.4f %10.4f %10.4f %10.4f\n",
+ vtz.vector(0),
+ vtz.vector(1),
+ r,
+ s.zc*0.1,
+ sum );
+ */
+ }
+ }
+ }
+ aida.cloud2D("HitGenNoScat/Z Hits vs cz",-1).fill( rktrk.cz(), nzmeas );
+ aida.cloud2D("HitGenNoScat/Nmeas vs cz",-1).fill( rktrk.cz(), nmeas );
+ aida.cloud2D("HitGenNoScat/Nz vs Nc",-1).fill( ncmeas, nzmeas );
+
+ return rkl;
+ }
+
+ // Generate perfect hits on the outward going arc only.
+ // Do not apply energy loss or multiple scattering.
+ public RKHitList GenerateMixedOneArcHits( RKTrack rktrk ){
+
+ // Convert to TRF style track.
+ toTRF trftrk = new toTRF(rktrk);
+ VTrack vtdca = trftrk.atDcaZaxis();
+
+ // This code only looks at the first arc of a track.
+ double pathlimit = trftrk.halfarc();
+
+ int nmeas = 0;
+ int nzmeas = 0;
+ int ncmeas = 0;
+ double sum = 0.;
+ RKHitList rkl = new RKHitList();
+ Propagator prop = rkgeom.newPropagator();
+
+ List<RKSurf> cyls = rkgeom.getCylinders();
+ List<RKSurf> zps = rkgeom.getZ( rktrk.z0(), rktrk.cz() );
+
+ int next_cyl = 0;
+ int next_zp = 0;
+ int max_cyl = cyls.size();
+ int max_zp = zps.size();
+
+ double sumde = 0;
+
+ // At top of the loop, vtdca contains the track parameters at the
+ // departure point for this step.
+ while ( (next_cyl<max_cyl) || ( next_zp < max_zp ) ){
+
+ VTrack vt_cyl = null;
+ VTrack vt_zp = null;
+
+ double path_cyl = 0.;
+ double path_zp = 0.;
+
+ int status_cyl = 0;
+ int status_zp = 0;
+
+ RKSurf rks_cyl = null;
+ RKSurf rks_zp = null;
+
+ PropStat ps_cyl = null;
+ PropStat ps_zp = null;
+
+ double r_cyl = 0.;
+ double z_cyl = 0.;
+ double r_zp = 0.;
+ double z_zp = 0.;
+
+ int next_cyl_sav = next_cyl;
+ int next_zp_sav = next_zp;
+
+ // Check cylindrical surfaces until the next good one is found.
+ //while ( next_cyl < -1 ){
+ while ( next_cyl < max_cyl ){
+
+ // Default status is failure.
+ status_cyl = 0;
+
+ // Parameters at the starting surface for this step.
+ vt_cyl = new VTrack(vtdca);
+
+ // Move to the next surface.
+ rks_cyl = cyls.get(next_cyl);
+ ps_cyl = prop.vecDirProp( vt_cyl, rks_cyl.getCylinder(), PropDir.FORWARD);
+
+ // If we do not get to this cylinder, try the next one.
+ // This can happen if starting point is outside this one but inside a later one.
+ if ( !ps_cyl.success() ) {
+ ++next_cyl;
+ continue;
+ }
+
+ // If we moved backwards, try the next cylinder outward.
+ path_cyl = ps_cyl.pathDistance();
+ if ( path_cyl < 0. ) {
+ //System.out.println ("Negative path at cyl..." );
+ status_cyl = 3;
+ ++next_cyl;
+ continue;
+ }
+ status_cyl = 1;
+
+ // Needed for downstream diagnostics.
+ r_cyl = rks_cyl.radius*mmTocm;
+ z_cyl = vt_cyl.vector(1);
+
+ if ( rks_cyl.inBounds(z_cyl) ){
+
+ // We have the next good surface from this list.
+ status_cyl = 2;
+ break;
+
+ } else{
+
+ // If not in bounds, try the next surface.
+ ++next_cyl;
+ continue;
+ }
+
+ }
+
+ // Check the next z surface.
+ while ( next_zp < max_zp ){
+
+ // Default status is failure.
+ status_zp = 0;
+
+ // Parameters at the start of this step.
+ vt_zp = new VTrack(vtdca);
+
+ // Move to the next surface.
+ rks_zp = zps.get(next_zp);
+ ps_zp = prop.vecDirProp( vt_zp, rks_zp.getZPlane(), PropDir.FORWARD );
+
+
+ // This happens when the z surface being tested is behind the starting point.
+ // Try the next surface.
+ if ( !ps_zp.success() ) {
+ ++next_zp;
+ continue;
+ }
+
+ // Normally this will not happen. If it does, try the next surface.
+ path_zp = ps_zp.pathDistance();
+ if ( path_zp < 0. ) {
+ status_zp = 3;
+ ++next_zp;
+ continue;
+ }
+ status_zp = 1;
+
+ // Needed for downstream diagnostics
+ r_zp = Math.sqrt(vt_zp.vector(0)*vt_zp.vector(0) +
+ vt_zp.vector(1)*vt_zp.vector(1) );
+ z_zp = rks_zp.zc*mmTocm;
+
+ if ( rks_zp.inBounds(r_zp) ){
+
+ // We have the next good surface from this list.
+ status_zp = 2;
+ break;
+ }else{
+
+ // Out of bounds so try the next surface.
+ ++next_zp;
+ continue;
+ }
+
+ }
+
+ // At this point we have zero, one or two good solutions.
+ // If two, choose the one with the shortest step length.
+ boolean addcyl = false;
+ boolean addzp = false;
+ if ( status_cyl == 2 && status_zp == 2 ){
+ if ( path_cyl < path_zp ) {
+ addcyl = true;
+ } else {
+ addzp = true;
+ }
+ } else if ( status_cyl == 2 ){
+ addcyl = true;
+
+ } else if ( status_zp == 2 ){
+ addzp = true;
+ } else {
+
+ // This should only happen when both lists are exhausted.
+ // Or when one of the hit types is turned off.
+ if ( next_cyl<max_cyl || next_zp < max_zp ){
+ System.out.println ("Don't know how I got here ... "
+ + next_cyl + " "
+ + max_cyl + " "
+ + next_zp + " "
+ + max_zp
+ );
+ }
+ }
+
+
+ // Add the chosen hit.
+ if ( addcyl ){
+
+ // Advance counter to next surface.
+ ++next_cyl;
+
+ // Restore the z index to it's position at the start of this step.
+ next_zp = next_zp_sav;
+
+ // Did we exceed the half-arc limit?
+ sum += path_cyl;
+ if ( sum > pathlimit ) break;
+
+ // Track parameters at this step.
+ vtdca = new VTrack(vt_cyl);
+
+ // Some diagnostics
+ aida.cloud2D( "HitGen/C r vs z",-1).fill( z_cyl, r_cyl );
+ aida.cloud2D( "HitGen/Both r vs z",-1).fill( z_cyl, r_cyl);
+
+ // Bookkeeping.
+ nmeas += rks_cyl.rodim;
+ ncmeas += rks_cyl.rodim;
+
+ if ( doms ){
+
+ // Interact.
+ // Should modify this to not interact if this is the last point on the outward trace.
+ VTrack test = new VTrack(vtdca);
+ if ( r_cyl < 1.3 ){
+ scatSim0.interact(vtdca);
+ } else if ( r_cyl < 10. ){
+ scatSim1.interact(vtdca);
+ }else{
+ scatSim2.interact(vtdca);
+ }
+
+
+ // Diagnostics.
+ double scatThick = 0.;
+ if ( r_cyl < 1.3 ){
+ scatThick = scatThick0;
+ }else if ( r_cyl < 10. ){
+ scatThick = scatThick1;
+ } else{
+ scatThick = scatThick2;
+ }
+ ThinCylMs scat1 = new ThinCylMs(scatThick);
+ TrackError vtmp = new TrackError();
+ ETrack et = new ETrack( test.surface().newPureSurface(), test.vector(), vtmp);
+ ETrack ets = new ETrack(et);
+ scat1.interact(et);
+
+ double dd = vtdca.vector().get(2)-test.vector().get(2);
+ double ee = et.error().get(2,2);
+ double r = dd/Math.sqrt(ee);
+ aida.cloud1D( "HitGen/Scat pull").fill(r);
+ aida.histogram1D("HitGen/Gen scat radius", 300, 0., 150.).fill(r_cyl);
+
+ }
+
+ // Add the hit.
+ rkl.addHit( new RKHit( rktrk, rks_cyl, vtdca, ps_cyl, sum ) );
+
+
+ } else if ( addzp ) {
+
+ // Advance counter to next surface.
+ ++next_zp;
+
+ // Restore the cylinder index to it's position at the start of this step.
+ next_cyl = next_cyl_sav;
+
+ // Did we exceed the half-arc limit?
+ sum += path_zp;
+ if ( sum > pathlimit ) break;
+
+ // Track parameters at this step.
+ vtdca = vt_zp;
+
+
+ // Some diagnostics
+ aida.cloud2D( "HitGen/Z r vs z",-1).fill( z_zp, r_zp );
+ aida.cloud2D( "HitGen/Both r vs z",-1).fill( z_zp, r_zp );
+
+ // Bookkeeping.
+ nmeas += rks_zp.rodim;
+ nzmeas += rks_zp.rodim;
+
+ if ( doms ){
+
+ // Interact.
+ // Should modify this to not interact if this is the last point on the outward trace.
+ VTrack test = new VTrack(vtdca);
+ if ( Math.abs(z_zp) < 25. ){
+ zscatSim1.interact(vtdca);
+ }else{
+ zscatSim2.interact(vtdca);
+ }
+ aida.histogram1D("HitGen/Gen scat z", 340, -170., 170.).fill(z_zp);
+
+ }
+
+ // Add the hit.
+ rkl.addHit( new RKHit( rktrk, rks_zp, vtdca, ps_zp, sum ) );
+
+
+ } else {
+
+ // Neither intersected but try next surfaces for both lists.
+ ++next_cyl;
+ ++next_zp;
+ }
+
+ }
+
+ // More diagnostics.
+ aida.cloud2D("HitGen/Path length vs cz",-1).fill(rktrk.cz(),sum);
+ aida.cloud2D("HitGen/C Hits vs cz",-1).fill( rktrk.cz(), ncmeas );
+ aida.cloud2D("HitGen/Z Hits vs cz",-1).fill( rktrk.cz(), nzmeas );
+ aida.cloud2D("HitGen/Nmeas vs cz",-1).fill( rktrk.cz(), nmeas );
+ aida.cloud2D("HitGen/Nz vs Nc",-1).fill( ncmeas, nzmeas );
+
+ aida.cloud1D("HitGen/Generated sumde:",-1).fill(sumde);
+ aida.histogram1D("HitGen/Generated sumde hist", 100, 0., 20.).fill(sumde*1000.);
+ RKDebug.Instance().setDeGen(sumde);
+
+ // If requested, also generate hits using the old algorithm and check that they
+ // agree with those from the new algorithm. This only makes sense if ms and eloss
+ // are disabled.
+ if ( doHitCheck && !( doeloss || doms ) ){
+ CheckHitList( rkl, rktrk );
+ }
+
+ return rkl;
+ }
+
+
+ // Check that the Mixed one arc hit generator matches the simple one.
+ // Only valid when scattering and eloss are turned off.
+ void CheckHitList( RKHitList hlist, RKTrack rktrk ){
+
+
+ RKHitList hlist2 = GenerateOneArcHits( rktrk );
+ int d1 = hlist.nHits() - hlist2.nHits();
+ int d2 = hlist.nCyl() - hlist2.nCyl();
+ int d3 = hlist.nZp() - hlist2.nZp();
+ int d4 = hlist.nDof() - hlist2.nDof();
+
+ if ( (d1==0) && (d2==0) && (d3==0) && (d4 ==0) ){
+ List<RKHit> l1 = hlist.getForward();
+ List<RKHit> l2 = hlist2.getForward();
+ for ( int i=0;i<l1.size(); ++i ){
+ RKHit h1 = l1.get(i);
+ RKHit h2 = l2.get(i);
+ double diff = h2.getPath()-h1.getPath();
+ aida.cloud1D("HitCheck/Difference in Path length: ").fill(diff);
+
+ }
+ }else{
+ System.out.println ( "Diffs: " + d1 + " " + d2 + " " + d3 + " " + d4 + " " + rktrk.cz() );
+ System.out.println ( "Hits1: "
+ + hlist2.nHits() + " "
+ + hlist2.nCyl() + " "
+ + hlist2.nZp() + " "
+ + hlist2.nDof() + " "
+ );
+ System.out.println ( "Hits2: "
+ + hlist.nHits() + " "
+ + hlist.nCyl() + " "
+ + hlist.nZp() + " "
+ + hlist.nDof() + " "
+ );
+ int i=0;
+ for ( RKHit h: hlist.getForward() ){
+ System.out.printf ( "Hits1: %3d %-20s %10.4f %10.4f\n",
+ i++,
+ h.getSurface().getTypeAsString(),
+ h.getSurface().radius,
+ h.getSurface().zc
+ );
+ }
+ i=0;
+ for ( RKHit h: hlist2.getForward() ){
+ System.out.printf ( "Hits2: %3d %-20s %10.4f %10.4f\n",
+ i++,
+ h.getSurface().getTypeAsString(),
+ h.getSurface().radius,
+ h.getSurface().zc
+ );
+ }
+
+ }
+
+ }
+
+}
+
