hps-java/src/main/java/org/lcsim/hps/recon/tracking
diff -u -r1.1 -r1.2
--- DataTrackerFakeHitDriver.java 9 Oct 2012 01:21:41 -0000 1.1
+++ DataTrackerFakeHitDriver.java 13 Oct 2012 02:57:30 -0000 1.2
@@ -21,7 +21,11 @@
import org.lcsim.event.SimTrackerHit;
import org.lcsim.event.base.BaseRawTrackerHit;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.fit.helicaltrack.HelixUtils;
+import org.lcsim.fit.helicaltrack.HitIdentifier;
import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.subdetector.BarrelEndcapFlag;
import org.lcsim.hps.users.phansson.WTrack;
import org.lcsim.hps.users.phansson.WTrackUtils;
import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHit;
@@ -47,7 +51,8 @@
WTrackUtils wutils = new WTrackUtils();
TrackerHitType trackerType = new TrackerHitType(TrackerHitType.CoordinateSystem.GLOBAL, TrackerHitType.MeasurementType.STRIP_1D);
CoordinateSystem coordinate_system = trackerType.getCoordinateSystem();
-
+ private HitIdentifier _ID = new HitIdentifier();
+ private boolean _doHth = false;
private String trackCollectionName = "MCParticle_HelicalTrackFit";
@@ -55,6 +60,7 @@
private String subdetectorName = "Tracker";
// Name of StripHit1D output collection.
private String stripHitOutputCollectionName = "StripClusterer_SiTrackerHitStrip1D";
+ private String hthOutputCollectionName = "RotatedHelicalTrackHits";
// Various data lists required by digitization.
private List<String> processPaths = new ArrayList<String>();
private List<IDetectorElement> processDEs = new ArrayList<IDetectorElement>();
@@ -67,12 +73,18 @@
private HashMap<SiSensor,IHistogram1D> _delta_itercount = new HashMap<SiSensor,IHistogram1D>();
IProfile1D _prf_final_deltas;
IProfile1D _prf_all_deltas;
+ IHistogram1D _h_nstriphits_top;
+ IHistogram1D _h_nstriphits_bottom;
IAnalysisFactory af = aida.analysisFactory();
IPlotter plotter_iter;
IPlotter plotter_itercount;
IPlotter plotter_iter_final;
IPlotter plotter_iter_all;
-
+
+ IPlotter plotter_nstripclusters;
+ IPlotter plotter_trackposodd;
+
+
int[][] counts = new int[2][10];
@@ -80,6 +92,10 @@
this.debug = debug;
}
+ public void setDoHth(boolean debug) {
+ this._doHth = debug;
+ }
+
// public void setReadoutCollectionName(String readoutCollectionName) {
// this.readoutCollectionName = readoutCollectionName;
// }
@@ -91,6 +107,10 @@
this.stripHitOutputCollectionName = stripHitOutputCollectionName;
}
+ public void setHthOutputCollectionName(String hthOutputCollectionName) {
+ this.hthOutputCollectionName = hthOutputCollectionName;
+ }
+
public void setTrackCollectionName(String trackCollectionName) {
this.trackCollectionName = trackCollectionName;
}
@@ -163,25 +183,30 @@
if(debug) System.out.println(this.getClass().getSimpleName() + ": found " + tracks.size() + " tracks (" + this.trackCollectionName + ")");
- if(tracks.isEmpty()) return;
+ //if(tracks.isEmpty()) return;
// Make new lists for output.
List<SiTrackerHit> stripHits1D = new ArrayList<SiTrackerHit>();
+ List<HelicalTrackHit> hths = new ArrayList<HelicalTrackHit>();
if(debug) System.out.println(this.getClass().getSimpleName() + ": Add hits for " + tracks.size() + " tracks (" + this.trackCollectionName + ")");
for(HelicalTrackFit helix : tracks) {
if(debug) System.out.println(this.getClass().getSimpleName() + ": trying to add hits for this track");
-
+
+ List<SiTrackerHit> stripHits1D_for_track = new ArrayList<SiTrackerHit>();
+ List<HelicalTrackHit> hths_for_track = new ArrayList<HelicalTrackHit>();
+
+
if(debug) {
System.out.println(this.getClass().getSimpleName() + helix.toString());
System.out.println(this.getClass().getSimpleName() + ": htf x0 " + helix.x0() + " y0 " + helix.y0());
+
+ System.out.println(this.getClass().getSimpleName() + ": create a WTrack object");
}
- if(debug) System.out.println(this.getClass().getSimpleName() + ": create a WTrack object");
-
- WTrack wtrack = new WTrack(helix,Math.abs(_bfield.z())); //remove sign from B-field (assumed to go along z-direction)
+ WTrack wtrack = new WTrack(helix,Math.abs(_bfield.z()),true); //remove sign from B-field (assumed to go along z-direction)
if(debug) {
System.out.println(this.getClass().getSimpleName() + ": " + wtrack.toString());
@@ -189,15 +214,22 @@
}
-
extrapolator.setTrack(helix.parameters());
+ int n_hits_top = 0;
+ int n_hits_bot = 0;
+ boolean isTopTrack = false;
+ boolean isBotTrack = false;
+
// Make hits if the helix passed through the sensor
for (SiSensor sensor : processSensors) {
- if(debug) System.out.println(this.getClass().getSimpleName() + ": add hits to sensor " + sensor.getName());
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": add hits to sensor " + sensor.getName() + " at position " + sensor.getGeometry().getPosition().toString());
- if(debug) System.out.println(this.getClass().getSimpleName() + ": sensor " + sensor.getName() + " position at " + sensor.getGeometry().getPosition().toString());
+ if(this._doHth && _ID.getLayer(sensor)%2==0) {
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": this was an even sensor -> skip for HTH production");
+ continue;
+ }
double zposition = sensor.getGeometry().getPosition().z();
@@ -214,23 +246,60 @@
if(isHit) {
if(debug) System.out.println(this.getClass().getSimpleName() + ": make a tracker hit and add to this sensor");
-
- stripHits1D.add(this.makeTrackerHit(sensor, wtrack));
+ if(SvtUtils.getInstance().isTopLayer(sensor)) {
+ n_hits_top++;
+ isTopTrack=true;
+ }
+ else {
+ n_hits_bot++;
+ isBotTrack=true;
+ }
+ if(this._doHth) {
+ hths_for_track.add(this.makeHelicalTrackHit(sensor, wtrack));
+ } else {
+ stripHits1D_for_track.add(this.makeTrackerHit(sensor, wtrack));
+
+ }
} else {
if(debug) System.out.println(this.getClass().getSimpleName() + ": this helix didn't pass within the sensor so no hit was added");
}
+
+
}
+
+ if(isTopTrack) this._h_nstriphits_top.fill(n_hits_top);
+ if(isBotTrack) this._h_nstriphits_bottom.fill(n_hits_bot);
+ if(isTopTrack && isBotTrack) {
+ System.out.println(this.getClass().getSimpleName() + ": tris track is both top and bottom??? \n" + wtrack.toString() + "\nHTF:" + helix.toString());
+ System.exit(1);
+ }
+
+ stripHits1D.addAll(stripHits1D_for_track);
+ hths.addAll(hths_for_track);
}
if (debug) {
- System.out.println("TrackerHit collection " + this.stripHitOutputCollectionName + " has " + stripHits1D.size() + " hits.");
+ System.out.println(this.getClass().getSimpleName() + ": Produced " + stripHits1D.size() + " hits for collection " + this.stripHitOutputCollectionName);
+ System.out.println(this.getClass().getSimpleName() + ": Produced " + hths.size() + " hits for collection " + this.hthOutputCollectionName);
}
// Put output hits into collection.
int flag = LCIOUtil.bitSet(0, 31, true); // Turn on 64-bit cell ID.
// event.put(this.rawTrackerHitOutputCollectionName, rawHits, RawTrackerHit.class, flag, toString());
event.put(this.stripHitOutputCollectionName, stripHits1D, SiTrackerHitStrip1D.class, 0, toString());
+ event.put("RotatedHelicalTrackHits", hths, HelicalTrackHit.class, 0);
+ //event.put("RotatedHelicalTrackHits", rotatedhits, HelicalTrackHit.class, 0);
+ if (debug) {
+ if(event.hasCollection(HelicalTrackHit.class, "RotatedHelicalTrackHits")) {
+ System.out.println(this.getClass().getSimpleName() + ": has hths:");
+ for(HelicalTrackHit hth : hths) {
+ System.out.println(this.getClass().getSimpleName() + ": " + hth.getPosition().toString());
+ }
+ } else {
+ System.out.println(this.getClass().getSimpleName() + ": has not hths!");
+ }
+ }
if (debug) {
for (int mod = 0; mod < 2; mod++) {
for (int layer = 0; layer < 10; layer++) {
@@ -307,17 +376,44 @@
//Rotate into tracking frame
Hep3Vector eta = VecOp.mult(detToTrk, w);
-
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": eta " + eta.toString());
+
Hep3Vector position = wutils.getHelixAndPlaneIntercept(wtrack, _orgloc, eta, h);
if(debug) System.out.println(this.getClass().getSimpleName() + ": found interception point at position " + position.toString());
+
+ HelicalTrackFit htf = wtrack._htf;
+ List<Double> s = HelixUtils.PathToXPlane(htf, position.x(), 0, 0);
+ Hep3Vector posOnHelix = HelixUtils.PointOnHelix(htf, s.get(0));
+ Hep3Vector posdiff = VecOp.sub(position, posOnHelix);
+ System.out.println(this.getClass().getSimpleName() + ": diffpos " + posdiff.toString() + " L " + position.toString() + " posOnHelix " + posOnHelix.toString() + " R=" + htf.R());
+
position = VecOp.mult(VecOp.inverse(detToTrk), position);
if(debug) System.out.println(this.getClass().getSimpleName() + ": rotate the hit position to the global frame -> " + position.toString());
+ // Need to make sure that the position is at the edge of the strip in the global frame
+ // 1. Rotate to the local sensor frame
+ position = ((SiSensor)electrodes.getDetectorElement()).getGeometry().getGlobalToLocal().transformed(position);
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": local (sensor) hit position " + position.toString());
+
+ // 2. Remove the coordinate of the unmeasured direction
+
+ position = new BasicHep3Vector(position.x(),0,position.z());
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": fixed local (sensor) hit position " + position.toString());
+
+ // 3. Transform back to global coordinates
+
+ position = ((SiSensor) electrodes.getDetectorElement()).getGeometry().getLocalToGlobal().transformed(position);
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": fixed global hit position " + position.toString());
+
+
//Fill some debug stuff
List<Double> deltas = wutils.getDeltas();
this._delta_itercount.get(sensor).fill(deltas.size());
@@ -328,13 +424,141 @@
}
_prf_final_deltas.fill(deltas.size(), deltas.get(deltas.size()-1));
-
+ //Fill dummy versions
SymmetricMatrix covariance = this.getCovariance(rth_cluster, electrodes);
- double time = getTime(rth_cluster);
- double energy = getEnergy(rth_cluster);
+ double time = this.getTime(rth_cluster);
+ double energy = this.getEnergy(rth_cluster);
SiTrackerHitStrip1D hit = new SiTrackerHitStrip1D(position, covariance, energy, time, rth_cluster, trackerType);
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": created SiStrip1D at " + position.toString());
+
+
+ return hit;
+
+
+ }
+
+
+
+
+ private HelicalTrackHit makeHelicalTrackHit(SiSensor sensor, WTrack wtrack) {
+ //private SiTrackerHitStrip1D makeTrackerHit(List<HPSFittedRawTrackerHit> cluster, SiSensorElectrodes electrodes) {
+ //create fake raw tracker hit
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": makeTrackerHit");
+
+ List<RawTrackerHit> rth_cluster = this.makeRawTrackerFakeHit(sensor);
+ if(rth_cluster.size()!=1) {
+ System.out.println(this.getClass().getSimpleName() + ": the fake raw tracker hit cluster is different than one!? " + rth_cluster.size() );
+ System.exit(1);
+ }
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": created a fake raw tracker hit ");
+ RawTrackerHit raw_hit = rth_cluster.get(0);
+ IIdentifier id = raw_hit.getIdentifier();
+
+ //Get the electrode objects
+ SiTrackerIdentifierHelper _sid_helper = (SiTrackerIdentifierHelper) sensor.getIdentifierHelper();
+ ChargeCarrier carrier = ChargeCarrier.getCarrier(_sid_helper.getSideValue(id));
+ SiSensorElectrodes electrodes = ((SiSensor) raw_hit.getDetectorElement()).getReadoutElectrodes(carrier);
+
+
+
+
+ ITransform3D local_to_global;
+ if (coordinate_system == TrackerHitType.CoordinateSystem.GLOBAL) {
+ local_to_global = new Transform3D();
+ }
+ else if (coordinate_system == TrackerHitType.CoordinateSystem.SENSOR) {
+ local_to_global = sensor.getGeometry().getLocalToGlobal();
+ }
+ else {
+ throw new RuntimeException(this.getClass().getSimpleName() + " problem with coord system " + coordinate_system.toString());
+ }
+
+
+ ITransform3D electrodes_to_global = electrodes.getLocalToGlobal();
+ ITransform3D global_to_hit = local_to_global.inverse();
+ ITransform3D electrodes_to_hit = Transform3D.multiply(global_to_hit,electrodes_to_global);
+
+ Hep3Vector u= electrodes_to_hit.rotated(electrodes.getMeasuredCoordinate(0));
+ Hep3Vector v= electrodes_to_hit.rotated(electrodes.getUnmeasuredCoordinate(0));
+ Hep3Vector w = VecOp.cross(u, v);
+ Hep3Vector _orgloc = new BasicHep3Vector(0., 0., 0.);
+ electrodes_to_global.transformed(_orgloc);
+ if(debug) {
+ System.out.println(this.getClass().getSimpleName() + ": electrodes u " + u.toString());
+ System.out.println(this.getClass().getSimpleName() + ": electrodes v " + v.toString());
+ System.out.println(this.getClass().getSimpleName() + ": electrodes w " + w.toString() + "( " + w.magnitude() + ")");
+ }
+
+ electrodes_to_global.getTranslation().translate(_orgloc);
+ if(debug) System.out.print(this.getClass().getSimpleName() + ": orgloc " + _orgloc.toString() + " -> ");
+ _orgloc = VecOp.mult(detToTrk, _orgloc);
+ if(debug) System.out.println(_orgloc.toString());
+
+
+
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": Try to find the interception point");
+
+ //B-field must go along Z direction
+ Hep3Vector h = new BasicHep3Vector(_bfield.x(),_bfield.y(),Math.abs(_bfield.z()));
+ h = VecOp.unit(h);
+ //Rotate into tracking frame
+ Hep3Vector eta = VecOp.mult(detToTrk, w);
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": eta " + eta.toString());
+
+ Hep3Vector position = wutils.getHelixAndPlaneIntercept(wtrack, _orgloc, eta, h);
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": found interception point at position " + position.toString());
+
+ if(debug) {
+ HelicalTrackFit htf = wtrack._htf;
+ List<Double> s = HelixUtils.PathToXPlane(htf, position.x(), 0, 0);
+ Hep3Vector posOnHelix = HelixUtils.PointOnHelix(htf, s.get(0));
+ Hep3Vector posdiff = VecOp.sub(position, posOnHelix);
+ System.out.println(this.getClass().getSimpleName() + ": diffpos " + posdiff.toString() + " L " + position.toString() + " posOnHelix " + posOnHelix.toString() + " R=" + htf.R());
+ }
+
+ //Fill some debug stuff
+ List<Double> deltas = wutils.getDeltas();
+ this._delta_itercount.get(sensor).fill(deltas.size());
+ IProfile1D dhisto = this._delta_histos.get(sensor);
+ for(int i=0;i<deltas.size();++i) {
+ this._prf_all_deltas.fill(i, deltas.get(i));
+ dhisto.fill(i,deltas.get(i));
+ }
+ _prf_final_deltas.fill(deltas.size(), deltas.get(deltas.size()-1));
+
+ //Fill dummy versions
+ //SymmetricMatrix covariance = this.getCovariance(rth_cluster, electrodes);
+ SymmetricMatrix covariance = new SymmetricMatrix(3);
+ covariance.setElement(0, 0, 1);
+ covariance.setElement(1, 1, 1);
+ covariance.setElement(2, 2, 1);
+
+ double time = this.getTime(rth_cluster);
+ double energy = this.getEnergy(rth_cluster);
+
+ //IDetectorElement de = sensor;
+ String det = _ID.getName(sensor);
+ int layer = _ID.getLayer(sensor);
+ BarrelEndcapFlag beflag = _ID.getBarrelEndcapFlag(sensor);
+
+ if(layer%2==0) {
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": problem, trying to create a HTH for even layer? " + layer);
+ System.exit(1);
+ }
+
+ HelicalTrackHit hit = new HelicalTrackHit(position, covariance,0.0, time, 3,rth_cluster, sensor.getName(), layer, beflag);
+ //SiTrackerHitStrip1D hit = new SiTrackerHitStrip1D(position, covariance, energy, time, rth_cluster, trackerType);
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": created HelicalTrackHit at " + position.toString());
+
+
return hit;
@@ -477,7 +701,12 @@
_prf_all_deltas = aida.profile1D("alldeltas",10,0,10);//,50,-20,20);
_prf_final_deltas = aida.profile1D("finaldeltas",10,0,10);//,50,-20,20);
-
+
+ _h_nstriphits_top = aida.histogram1D("NstripClusters top", 11, -0.5, 10.5);
+ _h_nstriphits_bottom = aida.histogram1D("NstripClusters bottom", 11, -0.5, 10.5);
+
+ //_h_trkposodd_top = aida.histogram1D("_h_trkposodd_top", 50, -0.5, 10.5);
+
plotter_iter_final = af.createPlotterFactory().create();
plotter_iter_final.createRegions(1,1);
plotter_iter_final.region(0).plot(_prf_final_deltas);
@@ -497,6 +726,16 @@
plotter_iter.style().statisticsBoxStyle().setVisible(false);
plotter_itercount = af.createPlotterFactory().create();
plotter_itercount.createRegions(5,4);
+
+ plotter_nstripclusters = af.createPlotterFactory().create();
+ plotter_nstripclusters.createRegions(2, 2);
+ plotter_nstripclusters.region(0).plot(_h_nstriphits_top);
+ plotter_nstripclusters.region(1).plot(_h_nstriphits_bottom);
+
+ plotter_trackposodd = af.createPlotterFactory().create();
+ plotter_trackposodd.createRegions(2, 2);
+ //plotter_trackposodd.region(0).plot(_h_trkposodd_top);
+
int i=0;
for(SiSensor sensor: this.processSensors) {
if(debug) System.out.println(this.getClass().getSimpleName() + ": " + i + ": adding plot to plotter for " + sensor.getName());
@@ -514,6 +753,7 @@
plotter_itercount.show();
plotter_iter_final.show();
plotter_iter_all.show();
+ plotter_nstripclusters.show();
}
}