hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl
diff -u -r1.7 -r1.8
--- GBLOutput.java 25 Aug 2013 21:42:19 -0000 1.7
+++ GBLOutput.java 26 Aug 2013 21:40:03 -0000 1.8
@@ -46,7 +46,8 @@
private MaterialSupervisor _materialmanager;
private MultipleScattering _scattering;
private Detector detector;
- private HPSTransformations _hpstrans = new HPSTransformations();
+ private HPSTransformations _hpstrans = new HPSTransformations();
+ private double _beamEnergy = 2.2; //GeV
@@ -123,9 +124,34 @@
// Use the sim tracker hits to find the truth particle
// this would give me the particle kinematics responsible for the track
+ List<MCParticle> mcp_pair = getAprimeDecayProducts(mcParticles);
+ if(_debug>0) {
+ double invMassTruth = Math.sqrt( Math.pow(mcp_pair.get(0).getEnergy()+mcp_pair.get(1).getEnergy(),2) - VecOp.add(mcp_pair.get(0).getMomentum(), mcp_pair.get(1).getMomentum()).magnitudeSquared() );
+ double invMassTruthTrks = getInvMassTracks(getHTF(mcp_pair.get(0)),getHTF(mcp_pair.get(1)));
+ System.out.printf("%s: invM = %f\n",this.getClass().getSimpleName(),invMassTruth);
+ System.out.printf("%s: invMTracks = %f\n",this.getClass().getSimpleName(),invMassTruthTrks);
+ }
MCParticle mcp = getMatchedTruthParticle(trk);
-
+
+ // cross-check
+ if(!mcp_pair.contains(mcp)) {
+ boolean hasBeamElectronParent = false;
+ for(MCParticle parent : mcp.getParents()) {
+ if(parent.getGeneratorStatus()!=MCParticle.FINAL_STATE && parent.getPDGID()==11 && parent.getMomentum().y()==0.0 && Math.abs(parent.getMomentum().magnitude() - _beamEnergy)<0.01) {
+ hasBeamElectronParent = true;
+ }
+ }
+ if(!hasBeamElectronParent) {
+ System.out.printf("%s: the matched MC particle is not an A' daughter and not a the recoil electrons!?\n",this.getClass().getSimpleName());
+ System.out.printf("%s: %s %d p %s org %s\n",this.getClass().getSimpleName(),mcp.getGeneratorStatus()==MCParticle.FINAL_STATE?"F":"I",mcp.getPDGID(),mcp.getMomentum().toString(),mcp.getOrigin().toString());
+ printMCParticles(mcParticles);
+ System.exit(1);
+ } else {
+ if(_debug>0) System.out.printf("%s: the matched MC particle is the recoil electron\n",this.getClass().getSimpleName());
+ }
+ }
+
List<SimTrackerHit> simHits = new ArrayList<SimTrackerHit>();
for (SimTrackerHit hit : simTrackerHits) {
if(hit.getMCParticle()==mcp) simHits.add(hit);
@@ -214,7 +240,7 @@
}
Hep3Vector trkposTruth = TrackUtils.getHelixPlaneIntercept(htfTruth, strip, Math.abs(_B.z()));
-
+
if(_debug>0) {
System.out.printf("trkposTruth at intercept %s\n",trkposTruth.toString());
}
@@ -247,7 +273,13 @@
System.out.printf("\n");
}
//System.exit(1);
- }
+ }
+ if(_debug>0) {
+ double s_truthSimHit = HelixUtils.PathToXPlane(htfTruth, simHit.getPositionVec().z(), 0, 0).get(0);
+ Hep3Vector trkposTruthSimHit = HelixUtils.PointOnHelix(htfTruth, s_truthSimHit);
+ Hep3Vector resTruthSimHit = VecOp.sub(this._hpstrans.transformVectorToTracking(simHit.getPositionVec()),trkposTruthSimHit);
+ System.out.printf("TruthSimHit residual %s for layer %d\n",resTruthSimHit.toString(),strip.layer());
+ }
//path length to intercept
double s = HelixUtils.PathToXPlane(htf,trkpos.x(),0,0).get(0);
@@ -649,9 +681,74 @@
return chi2.e(0, 0);
}
+
+
+ private List<MCParticle> getAprimeDecayProducts(List<MCParticle> mcParticles) {
+ List<MCParticle> pair = new ArrayList<MCParticle>();
+ for(MCParticle mcp : mcParticles) {
+ if(mcp.getGeneratorStatus()!=MCParticle.FINAL_STATE) continue;
+ boolean hasAprimeParent = false;
+ for(MCParticle parent : mcp.getParents()) {
+ if(Math.abs(parent.getPDGID())==622) hasAprimeParent = true;
+ }
+ if(hasAprimeParent) pair.add(mcp);
+ }
+ if(pair.size()!=2) {
+ System.out.printf("%s: ERROR this event has %d mcp with 622 as parent!!?? \n",this.getClass().getSimpleName(),pair.size());
+ this.printMCParticles(mcParticles);
+ System.exit(1);
+ }
+ if( Math.abs(pair.get(0).getPDGID()) != 11 || Math.abs(pair.get(1).getPDGID()) != 11 ) {
+ System.out.printf("%s: ERROR decay products are not e+e-? \n",this.getClass().getSimpleName());
+ this.printMCParticles(mcParticles);
+ System.exit(1);
+ }
+ if(pair.get(0).getPDGID()*pair.get(1).getPDGID() > 0) {
+ System.out.printf("%s: ERROR decay products have the same sign? \n",this.getClass().getSimpleName());
+ this.printMCParticles(mcParticles);
+ System.exit(1);
+ }
+ return pair;
+
+ }
+
+ private void printMCParticles(List<MCParticle> mcParticles) {
+ System.out.printf("%s: printMCParticles \n",this.getClass().getSimpleName());
+ System.out.printf("%s: %d mc particles \n",this.getClass().getSimpleName(),mcParticles.size());
+ for(MCParticle mcp : mcParticles) {
+ if(mcp.getGeneratorStatus()!=MCParticle.FINAL_STATE) continue;
+ System.out.printf("\n%s: (%s) %d p %s org %s %s \n",this.getClass().getSimpleName(),
+ mcp.getGeneratorStatus()==MCParticle.FINAL_STATE?"F":"I",mcp.getPDGID(),mcp.getMomentum().toString(),mcp.getOrigin().toString(),
+ mcp.getParents().size()>0?"parents:":"");
+ for(MCParticle parent : mcp.getParents()) {
+ System.out.printf("%s: (%s) %d p %s org %s %s \n",this.getClass().getSimpleName(),
+ parent.getGeneratorStatus()==MCParticle.FINAL_STATE?"F":"I",parent.getPDGID(),parent.getMomentum().toString(),parent.getOrigin().toString(),
+ parent.getParents().size()>0?"parents:":"");
+ for(MCParticle grparent : parent.getParents()) {
+ System.out.printf("%s: (%s) %d p %s org %s %s \n",this.getClass().getSimpleName(),
+ grparent.getGeneratorStatus()==MCParticle.FINAL_STATE?"F":"I",grparent.getPDGID(),grparent.getMomentum().toString(),grparent.getOrigin().toString(),
+ grparent.getParents().size()>0?"parents:":"");
+ }
+
+ }
+ }
+ return;
+ }
-
+ private double getInvMassTracks(HelicalTrackFit htf1, HelicalTrackFit htf2) {
+ double p1 = htf1.p(this._B.magnitude());
+ double p2 = htf2.p(this._B.magnitude());
+ Hep3Vector p1vec = VecOp.mult(p1, HelixUtils.Direction(htf1, 0));
+ Hep3Vector p2vec = VecOp.mult(p2, HelixUtils.Direction(htf2, 0));
+ double me = 0.000510998910;
+ double E1 = Math.sqrt(p1*p1 + me*me);
+ double E2 = Math.sqrt(p2*p2 + me*me);
+ //System.out.printf("p1 %f %s E1 %f\n",p1,p1vec.toString(),E1);
+ //System.out.printf("p2 %f %s E2 %f\n",p2,p2vec.toString(),E2);
+ return Math.sqrt( Math.pow(E1+E2,2) - VecOp.add(p1vec, p2vec).magnitudeSquared() );
+ }
+
public class PerigeeParams {
private BasicMatrix _params;