Author: [log in to unmask]
Date: Wed Nov 18 14:35:46 2015
New Revision: 3968
Log:
Fix MC truth processing.
Modified:
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutputDriver.java
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java Wed Nov 18 14:35:46 2015
@@ -58,7 +58,6 @@
private final MaterialSupervisor materialManager;
private final MultipleScattering _scattering;
private final double _beamEnergy = 1.1; //GeV
- private boolean AprimeEvent = false; // do extra checks
private boolean hasXPlanes = false;
private boolean addBeamspot = false;
private double beamspotTiltZOverY = 0; //Math.PI/180* 15;
@@ -140,10 +139,6 @@
}
}
- void setAPrimeEventFlag(boolean flag) {
- this.AprimeEvent = flag;
- }
-
void setXPlaneFlag(boolean flag) {
this.hasXPlanes = flag;
}
@@ -170,25 +165,61 @@
// Find the truth particle of the track
MCParticle mcp = null;
-
+ MCParticle ap = null;
+
+ // MC processing
if (isMC) {
+
+ // find the truth particle for this track
mcp = getMatchedTruthParticle(trk);
+ // check if this is an A' event
+ for(MCParticle part : mcParticles) {
+ if(Math.abs(part.getPDGID()) == 622) {
+ ap = part;
+ break;
+ }
+ }
+
if (mcp == null) {
System.out.printf("%s: WARNING!! no truth particle found in event!\n", this.getClass().getSimpleName());
this.printMCParticles(mcParticles);
//System.exit(1);
- } else if (_debug > 0) {
- System.out.printf("%s: truth particle (pdgif %d ) found in event!\n", this.getClass().getSimpleName(), mcp.getPDGID());
- }
-
- if (AprimeEvent) {
- checkAprimeTruth(mcp, mcParticles);
- }
- }
+ } else {
+ if (_debug > 0) System.out.printf("%s: truth particle (pdgif %d ) found in event!\n", this.getClass().getSimpleName(), mcp.getPDGID());
+
+ // If this is an A' event, do some more checks
+ if ( ap != null) {
+ // A few MC files have broken links b/w parents-daughters
+ // This causes the MC particle to come from the origin even if the decay happen somewhere else
+ if(this.getAprimeDecayProducts(mcParticles).size()>0) {
+ //do a full check
+ checkAprimeTruth(mcp, mcParticles);
+ }
+ }
+ }
+ }
// Get track parameters from MC particle
- HelicalTrackFit htfTruth = (isMC && mcp != null) ? TrackUtils.getHTF(mcp, -1.0 * this.bFieldVector.z()) : null;
+ HelicalTrackFit htfTruth = null;
+
+ if( isMC && mcp != null) {
+ // check if we should be using a different origin than the particle tells us
+ Hep3Vector mcp_origin;
+ if( ap != null) {
+ // There is an A' here. Use its origin if different
+ if (_debug > 0) System.out.printf("%s: A' found with origin %s compared to particle %s (diff: %s)\n", this.getClass().getSimpleName(), ap.getOrigin().toString(), mcp.getOrigin().toString(), VecOp.sub(ap.getOrigin(), mcp.getOrigin()).toString());
+ if(VecOp.sub(ap.getOrigin(), mcp.getOrigin()).magnitude() > 0.00001)
+ mcp_origin = ap.getOrigin();
+ else
+ mcp_origin = mcp.getOrigin();
+ } else {
+ // No A', use particle origin
+ mcp_origin = mcp.getOrigin();
+ }
+
+ htfTruth = TrackUtils.getHTF(mcp,mcp_origin, -1.0 * this.bFieldVector.z());
+ }
// Use the truth helix as the initial track for GBL?
//htf = htfTruth;
@@ -677,11 +708,53 @@
}
+
+ 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);
+ }
+ }
+
+ return pair;
+
+ }
+
private void checkAprimeTruth(MCParticle mcp, List<MCParticle> mcParticles) {
+
List<MCParticle> mcp_pair = getAprimeDecayProducts(mcParticles);
+
+ if (mcp_pair.size() != 2) {
+ System.out.printf("%s: ERROR this event has %d mcp with 622 as parent!!?? \n", this.getClass().getSimpleName(), mcp_pair.size());
+ this.printMCParticles(mcParticles);
+ System.exit(1);
+ }
+ if (Math.abs(mcp_pair.get(0).getPDGID()) != 11 || Math.abs(mcp_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 (mcp_pair.get(0).getPDGID() * mcp_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);
+ }
+
+
+
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(TrackUtils.getHTF(mcp_pair.get(0), -1.0 * this.bFieldVector.z()), TrackUtils.getHTF(mcp_pair.get(1), -1.0 * this.bFieldVector.z()));
@@ -896,40 +969,6 @@
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());
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutputDriver.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutputDriver.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutputDriver.java Wed Nov 18 14:35:46 2015
@@ -79,7 +79,6 @@
gbl = new GBLOutput(gblFileName, bfield); // if filename is empty no text file is written
gbl.setDebug(_debug);
gbl.buildModel(detector);
- gbl.setAPrimeEventFlag(false);
gbl.setXPlaneFlag(false);
gbl.setAddBeamspot(addBeamspot);
gbl.setBeamspotScatAngle(beamspotScatAngle);
