1 package org.lcsim.util.heprep;
2
3 import static java.lang.Math.abs;
4 import hep.graphics.heprep.HepRepFactory;
5 import hep.graphics.heprep.HepRepInstance;
6 import hep.graphics.heprep.HepRepInstanceTree;
7 import hep.graphics.heprep.HepRepType;
8 import hep.graphics.heprep.HepRepTypeTree;
9 import hep.physics.particle.properties.UnknownParticleIDException;
10 import hep.physics.vec.BasicHep3Vector;
11 import hep.physics.vec.Hep3Vector;
12 import hep.physics.vec.VecOp;
13
14 import java.awt.Color;
15 import java.util.List;
16
17 import org.lcsim.event.EventHeader;
18 import org.lcsim.event.EventHeader.LCMetaData;
19 import org.lcsim.event.MCParticle;
20 import org.lcsim.geometry.Detector;
21 import org.lcsim.util.swim.HelixSwimmer;
22
23
24
25
26
27
28
29
30
31 class MCParticleConverter implements HepRepCollectionConverter
32 {
33 private static final double[] IP = { 0,0,0 };
34
35
36 public boolean canHandle(Class k)
37 {
38 return MCParticle.class.isAssignableFrom(k);
39 }
40 public void convert(EventHeader event, List collection, HepRepFactory factory, HepRepTypeTree typeTree, HepRepInstanceTree instanceTree)
41 {
42 LCMetaData meta = event.getMetaData(collection);
43 String name = meta.getName();
44
45 Detector detector = null;
46
47
48
49 try {
50 detector = event.getDetector();
51 }
52 catch (Exception x)
53 {}
54
55 double trackingRMax = 10000;
56
57
58 if (detector != null) {
59 trackingRMax = detector.getConstants().get("tracking_region_radius").getValue();
60
61 }
62
63 double ptMinCut = 0.05;
64 double rCut = 1.0;
65 double[] field;
66 if (detector != null)
67 field = detector.getFieldMap().getField(IP);
68 else
69 field = new double[3];
70
71 HelixSwimmer helix = new HelixSwimmer(field[2]);
72
73 HepRepType typeX = factory.createHepRepType(typeTree, name);
74 typeX.addAttValue("layer",LCSimHepRepConverter.PARTICLES_LAYER);
75 typeX.addAttValue("drawAs","Line");
76
77 typeX.addAttDef("momentum","Particle Momentum", "physics", "GeV");
78 typeX.addAttDef("energy","Particle Energy","physcs","GeV");
79 typeX.addAttDef("pT","Particle Transverse Energy","physics","GeV");
80 typeX.addAttDef("time","Particle Production Time","physics","nanoseconds");
81 typeX.addAttDef("type","Particle Type", "physics", "");
82
83
84 HepRepType neutralType = factory.createHepRepType(typeX, "Neutral");
85 neutralType.addAttValue("color",Color.ORANGE);
86
87 HepRepType photonType = factory.createHepRepType(neutralType, "Photon");
88 photonType.addAttValue("color",Color.YELLOW);
89
90 HepRepType neutrinoType = factory.createHepRepType(neutralType, "Neutrino");
91 neutrinoType.addAttValue("color",Color.ORANGE);
92
93 HepRepType neutralHadronType = factory.createHepRepType(neutralType, "Neutral hadron");
94 neutralHadronType.addAttValue("color",Color.GREEN);
95
96 HepRepType chargedType = factory.createHepRepType(typeX, "Charged");
97 chargedType.addAttValue("color",Color.BLUE);
98
99 HepRepInstance charged = factory.createHepRepInstance(instanceTree, typeX);
100 HepRepInstance neutral = factory.createHepRepInstance(instanceTree, typeX);
101
102 for (MCParticle p : (List<MCParticle>) collection)
103 {
104 try
105 {
106 Hep3Vector start = p.getOrigin();
107 Hep3Vector momentum = p.getMomentum();
108 double charge = p.getCharge();
109 helix.setTrack(momentum, start, (int) charge);
110 Hep3Vector stop;
111
112 try
113 {
114 stop = p.getEndPoint();
115
116 if (stop.x() == 0 && stop.y() == 0 && stop.z() == 0)
117 {
118 if(p.getGeneratorStatus()==MCParticle.FINAL_STATE) stop = helix.getPointAtDistance(trackingRMax);
119 }
120 }
121 catch (RuntimeException x)
122 {
123
124 if(p.getGeneratorStatus()==MCParticle.FINAL_STATE)
125 {
126 stop = helix.getPointAtDistance(trackingRMax);
127 }
128 else
129 {
130 stop = new BasicHep3Vector();
131 }
132 }
133
134 if (charge == 0 || field[2] == 0)
135 {
136 HepRepInstance instanceX = factory.createHepRepInstance(charged, chargedType);
137 if(charge == 0)
138 {
139 int pdgId = p.getPDGID();
140
141 HepRepType type = neutralHadronType;
142 if(isNeutrino(pdgId)) type = neutrinoType;
143 if(abs(pdgId)==22) type = photonType;
144 instanceX = factory.createHepRepInstance(neutral, type);
145 }
146
147 setDefaultAttValues(instanceX,p);
148
149 factory.createHepRepPoint(instanceX,start.x(),start.y(),start.z());
150 factory.createHepRepPoint(instanceX,stop.x(),stop.y(),stop.z());
151
152 }
153 else
154 {
155 double pT = Math.sqrt(momentum.x()*momentum.x()+momentum.y()*momentum.y());
156
157 double r = Math.sqrt(start.x()*start.x()+start.y()*start.y());
158 if(pT>ptMinCut || (pT<ptMinCut && r>rCut))
159 {
160 double dAlpha = 10;
161 HepRepInstance instanceX = factory.createHepRepInstance(charged, chargedType);
162
163 setDefaultAttValues(instanceX,p);
164
165
166 factory.createHepRepPoint(instanceX,start.x(),start.y(),start.z());
167 double absZ = Math.abs(stop.z());
168 double rSquared = stop.x()*stop.x()+stop.y()*stop.y();
169 Hep3Vector point = start;
170
171 for (int k = 1;k<200;k++)
172 {
173 double d = VecOp.sub(point,stop).magnitudeSquared();
174
175 if (d < 2)
176 {
177 factory.createHepRepPoint(instanceX,stop.x(),stop.y(),stop.z());
178 break;
179 }
180 else if (Math.abs(point.z()) > absZ ||
181 point.x()*point.x()+point.y()*point.y() > rSquared)
182 {
183 break;
184 }
185 else
186 {
187 point = helix.getPointAtDistance(k*dAlpha);
188 factory.createHepRepPoint(instanceX,point.x(),point.y(),point.z());
189 }
190 }
191 }
192 }
193 }
194 catch (UnknownParticleIDException x)
195 {
196
197 }
198 }
199
200
201
202
203
204
205
206
207
208
209
210
211 }
212
213 boolean isNeutrino(int pdgId)
214 {
215 if(abs(pdgId)==12) return true;
216 if(abs(pdgId)==14) return true;
217 if(abs(pdgId)==16) return true;
218
219 return false;
220 }
221
222
223 private void setDefaultAttValues(HepRepInstance instanceX, MCParticle p)
224 {
225 double x = p.getMomentum().x();
226 double y = p.getMomentum().y();
227 double pT = Math.sqrt(x*x + y*y);
228
229 instanceX.addAttValue("pT",pT);
230 instanceX.addAttValue("particle",p.getType().getName());
231 instanceX.addAttValue("energy",p.getEnergy());
232 instanceX.addAttValue("momentum",p.getMomentum().magnitude());
233 instanceX.addAttValue("time",p.getProductionTime());
234 }
235
236 }