Skip to content

Commit 439b49d

Browse files
author
Matthew Solt
committed
adding MC Full Detector Truth File
1 parent fbfb2c7 commit 439b49d

1 file changed

Lines changed: 321 additions & 0 deletions

File tree

‎analysis/src/main/java/org/hps/analysis/MC/MCFullDetectorTruth.java‎

Lines changed: 321 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,321 @@
1+
package org.hps.analysis.MC;
2+
3+
import hep.physics.vec.BasicHep3Vector;
4+
import hep.physics.vec.Hep3Vector;
5+
import hep.physics.vec.VecOp;
6+
7+
import java.util.ArrayList;
8+
import java.util.HashMap;
9+
import java.util.List;
10+
import java.util.Map;
11+
12+
import org.hps.recon.tracking.CoordinateTransformations;
13+
import org.hps.recon.tracking.TrackUtils;
14+
import org.lcsim.event.MCParticle;
15+
import org.lcsim.event.SimCalorimeterHit;
16+
import org.lcsim.event.SimTrackerHit;
17+
import org.lcsim.geometry.FieldMap;
18+
import org.lcsim.geometry.IDDecoder;
19+
import org.lcsim.util.swim.Trajectory;
20+
21+
/**
22+
* This is the tuple template driver
23+
* Use this to add your code and variables to make a tuple
24+
* Run the GeneralTupleDriver to output info into a text file
25+
* Change the steering file to include this driver
26+
* Run "makeTree.py" on text file to create a root tuple
27+
*
28+
* @author mrsolt on Aug 31, 2017
29+
*/
30+
31+
public class MCFullDetectorTruth{
32+
33+
public static Map<MCParticle, List<SimTrackerHit>> BuildTrackerHitMap(List<SimTrackerHit> trackerHits){
34+
Map<MCParticle, List<SimTrackerHit>> trackerHitMap = new HashMap<MCParticle, List<SimTrackerHit>>();
35+
for (SimTrackerHit hit : trackerHits) {
36+
MCParticle p = hit.getMCParticle();
37+
if (p == null) {
38+
throw new RuntimeException("Tracker hit points to null MCParticle!");
39+
}
40+
if (trackerHitMap.get(p) == null) {
41+
trackerHitMap.put(p, new ArrayList<SimTrackerHit>());
42+
}
43+
trackerHitMap.get(p).add(hit);
44+
}
45+
return trackerHitMap;
46+
}
47+
48+
public static Map<MCParticle, List<SimCalorimeterHit>> BuildCalHitMap(List<SimCalorimeterHit> calHits){
49+
// map particle to a list of its sim cal hits
50+
Map<MCParticle, List<SimCalorimeterHit>> calHitMap = new HashMap<MCParticle, List<SimCalorimeterHit>>();
51+
for (SimCalorimeterHit hit : calHits) {
52+
int nmc = hit.getMCParticleCount();
53+
for (int i = 0; i < nmc; i++) {
54+
MCParticle p = hit.getMCParticle(i);
55+
if (p == null) {
56+
throw new RuntimeException("Cal hit points to null MCParticle!");
57+
}
58+
if (calHitMap.get(p) == null) {
59+
calHitMap.put(p, new ArrayList<SimCalorimeterHit>());
60+
}
61+
if (!calHitMap.get(p).contains(hit)) {
62+
calHitMap.get(p).add(hit);
63+
}
64+
}
65+
}
66+
return calHitMap;
67+
}
68+
69+
public static String trackHitLayer(IDDecoder trackerDecoder,SimTrackerHit hit) {
70+
String layer = Integer.toString( (int) hit.getLayer());
71+
double y = hit.getPositionVec().y();
72+
String volume = "";
73+
if(y > 0) volume = "t";
74+
else volume = "b";
75+
String prefix = "L" + layer + volume;
76+
return prefix;
77+
}
78+
79+
public static String MCParticleType(MCParticle p, List<MCParticle> particles, double ebeam) {
80+
boolean isEle = false;
81+
boolean isPos = false;
82+
83+
// particle PDG ID
84+
int pdgid = p.getPDGID();
85+
86+
// particle energy
87+
double energy = p.getEnergy();
88+
89+
if(pdgid == 11){
90+
isEle = true;
91+
}
92+
93+
if(pdgid == -11){
94+
isPos = true;
95+
}
96+
97+
// parent index in particle list
98+
Integer parIdx = null;
99+
100+
// find parent's index in the particle collection
101+
if (p.getParents().size() > 0) {
102+
MCParticle parent = p.getParents().get(0);
103+
for (int i = 0; i < particles.size(); i++) {
104+
if (particles.get(i) == parent) {
105+
parIdx = i;
106+
break;
107+
}
108+
}
109+
}
110+
else{
111+
parIdx = -1;
112+
}
113+
if(isEle && parIdx == 0){
114+
for (MCParticle particle : particles) {
115+
if(particle.getParents().size() == 0) continue;
116+
if(particle.getPDGID() == 11 && particle.getParents().get(0).getPDGID() == 622 && p != particle){
117+
if(particle.getEnergy() < energy){
118+
return "triEle1";
119+
}
120+
else{
121+
return "triEle2";
122+
}
123+
}
124+
}
125+
}
126+
127+
if(isPos && parIdx == 0){
128+
return "triPos";
129+
}
130+
131+
132+
MCParticle wab = null;
133+
134+
MCParticle ele1 = null;// conversion electron daughter
135+
MCParticle ele2 = null;// recoil wab electron
136+
MCParticle pos = null;// conversion positron daughter
137+
138+
List<MCParticle> wabParticles = null;
139+
140+
for (MCParticle particle : particles) {
141+
if (particle.getPDGID() == 22) {
142+
if (particle.getDaughters().size() != 2) continue;
143+
double wabEnergy = particle.getEnergy();
144+
for(MCParticle part : particles){
145+
if(part.getPDGID() != 11 || !part.getParents().isEmpty()) continue;
146+
double eleEnergy = part.getEnergy();
147+
double esum = wabEnergy + eleEnergy;
148+
if(esum < 0.98 * ebeam || esum > 1.02 * ebeam) continue;
149+
ele2 = part;
150+
wab = particle;
151+
wabParticles = wab.getDaughters();
152+
break;
153+
}
154+
}
155+
}
156+
if (wab == null) {
157+
return "";
158+
}
159+
160+
161+
for (MCParticle particle : wabParticles) {
162+
if(particle.getPDGID() == 11){
163+
pos = particle;
164+
}
165+
if(particle.getPDGID() == -11){
166+
ele1 = particle;
167+
}
168+
}
169+
170+
if (ele1 == p) {
171+
return "wabEle1";
172+
}
173+
if (ele2 == p) {
174+
return "wabEle2";
175+
}
176+
if (pos == p) {
177+
return "wabPos";
178+
}
179+
180+
return "";
181+
}
182+
183+
public static Hep3Vector extrapolateTrackPosition(Hep3Vector startPosition, Hep3Vector startMomentum, double endPositionZ, double stepSize, FieldMap fieldMap, double q) {
184+
185+
// Start by transforming detector vectors into tracking frame
186+
Hep3Vector currentPosition = CoordinateTransformations.transformVectorToTracking(startPosition);
187+
Hep3Vector currentMomentum = CoordinateTransformations.transformVectorToTracking(startMomentum);
188+
189+
//System.out.println(currentPosition.x() + " " + currentPosition.y() + " " + currentPosition.z() + " ");
190+
191+
// Retrieve the y component of the bfield in the middle of detector
192+
double bFieldY = fieldMap.getField(new BasicHep3Vector(0, 0, 500.0)).y();
193+
194+
// HACK: LCSim doesn't deal well with negative fields so they are
195+
// turned to positive for tracking purposes. As a result,
196+
// the charge calculated using the B-field, will be wrong
197+
// when the field is negative and needs to be flipped.
198+
if (bFieldY < 0)
199+
q = q * (-1);
200+
201+
// Swim the track through the B-field until the end point is reached.
202+
// The position of the track will be incremented according to the step
203+
// size up to ~90% of the final position. At this point, a finer
204+
// track size will be used.
205+
boolean stepSizeChange = false;
206+
while (currentPosition.x() < endPositionZ) {
207+
// The field map coordinates are in the detector frame so the
208+
// extrapolated track position needs to be transformed from the
209+
// track frame to detector.
210+
Hep3Vector currentPositionDet = CoordinateTransformations.transformVectorToDetector(currentPosition);
211+
212+
// Get the field at the current position along the track.
213+
bFieldY = fieldMap.getField(currentPositionDet).y();
214+
215+
// Get a trajectory (Helix or Line objects) created with the
216+
// track parameters at the current position.
217+
Trajectory trajectory = TrackUtils.getTrajectory(currentMomentum, new org.lcsim.spacegeom.SpacePoint(currentPosition), q, bFieldY);
218+
219+
// Using the new trajectory, extrapolated the track by a step and
220+
// update the extrapolated position.
221+
//if(CoordinateTransformations.transformVectorToTracking(startPosition).x() > currentPosition.x()){
222+
if(currentMomentum.x() < 0){
223+
//System.out.println("Extrapolation is going backwards! Break loop! Z = " + currentPosition.x());
224+
break;
225+
}
226+
currentPosition = trajectory.getPointAtDistance(stepSize);
227+
228+
// Calculate the momentum vector at the new position. This will
229+
// be used when creating the trajectory that will be used to
230+
// extrapolate the track in the next iteration.
231+
currentMomentum = VecOp.mult(currentMomentum.magnitude(), trajectory.getUnitTangentAtLength(stepSize));
232+
233+
//System.out.println(currentPosition.x() + " " + currentPosition.y() + " " + currentPosition.z() + " ");
234+
235+
// If the position of the track along X (or z in the detector frame)
236+
// is at 90% of the total distance, reduce the step size.
237+
if (currentPosition.x() / endPositionZ > .80 && !stepSizeChange) {
238+
stepSize /= 10;
239+
// System.out.println("Changing step size: " + stepSize);
240+
stepSizeChange = true;
241+
}
242+
}
243+
244+
// Transform vector back to detector frame
245+
return CoordinateTransformations.transformVectorToDetector(currentPosition);
246+
}
247+
248+
public static Hep3Vector extrapolateTrackMomentum(Hep3Vector startPosition, Hep3Vector startMomentum, double endPositionZ, double stepSize, FieldMap fieldMap, double q) {
249+
250+
// Start by transforming detector vectors into tracking frame
251+
Hep3Vector currentPosition = CoordinateTransformations.transformVectorToTracking(startPosition);
252+
Hep3Vector currentMomentum = CoordinateTransformations.transformVectorToTracking(startMomentum);
253+
254+
// Retrieve the y component of the bfield in the middle of detector
255+
double bFieldY = fieldMap.getField(new BasicHep3Vector(0, 0, 500.0)).y();
256+
257+
// HACK: LCSim doesn't deal well with negative fields so they are
258+
// turned to positive for tracking purposes. As a result,
259+
// the charge calculated using the B-field, will be wrong
260+
// when the field is negative and needs to be flipped.
261+
if (bFieldY < 0)
262+
q = q * (-1);
263+
264+
// Swim the track through the B-field until the end point is reached.
265+
// The position of the track will be incremented according to the step
266+
// size up to ~90% of the final position. At this point, a finer
267+
// track size will be used.
268+
boolean stepSizeChange = false;
269+
while (currentPosition.x() < endPositionZ) {
270+
// The field map coordinates are in the detector frame so the
271+
// extrapolated track position needs to be transformed from the
272+
// track frame to detector.
273+
Hep3Vector currentPositionDet = CoordinateTransformations.transformVectorToDetector(currentPosition);
274+
275+
// Get the field at the current position along the track.
276+
bFieldY = fieldMap.getField(currentPositionDet).y();
277+
278+
// Get a trajectory (Helix or Line objects) created with the
279+
// track parameters at the current position.
280+
Trajectory trajectory = TrackUtils.getTrajectory(currentMomentum, new org.lcsim.spacegeom.SpacePoint(currentPosition), q, bFieldY);
281+
282+
//if(CoordinateTransformations.transformVectorToTracking(startPosition).x() > currentPosition.x()){
283+
if(currentMomentum.x() < 0){
284+
//System.out.println("Extrapolation is going backwards! Break loop! Z = " + currentPosition.x());
285+
break;
286+
}
287+
288+
// Using the new trajectory, extrapolated the track by a step and
289+
// update the extrapolated position.
290+
currentPosition = trajectory.getPointAtDistance(stepSize);
291+
292+
// Calculate the momentum vector at the new position. This will
293+
// be used when creating the trajectory that will be used to
294+
// extrapolate the track in the next iteration.
295+
currentMomentum = VecOp.mult(currentMomentum.magnitude(), trajectory.getUnitTangentAtLength(stepSize));
296+
297+
// If the position of the track along X (or z in the detector frame)
298+
// is at 90% of the total distance, reduce the step size.
299+
if (currentPosition.x() / endPositionZ > .80 && !stepSizeChange) {
300+
stepSize /= 10;
301+
// System.out.println("Changing step size: " + stepSize);
302+
stepSizeChange = true;
303+
}
304+
}
305+
306+
// Transform vector back to detector frame
307+
return CoordinateTransformations.transformVectorToDetector(currentMomentum);
308+
}
309+
310+
public static double deltaThetaX(Hep3Vector p1, Hep3Vector p2){
311+
double theta1 = p1.x()/p1.z();
312+
double theta2 = p2.x()/p2.z();
313+
return theta2 - theta1;
314+
}
315+
316+
public static double deltaThetaY(Hep3Vector p1, Hep3Vector p2){
317+
double theta1 = p1.y()/p1.z();
318+
double theta2 = p2.y()/p2.z();
319+
return theta2 - theta1;
320+
}
321+
}

0 commit comments

Comments
 (0)