1 package org.lcsim.mc.fast.util;
2
3 import hep.physics.vec.Hep3Vector;
4 import java.util.ArrayList;
5 import java.util.List;
6 import org.lcsim.event.MCParticle;
7 import org.lcsim.util.Driver;
8 import org.lcsim.event.EventHeader;
9 import static org.lcsim.mc.fast.util.MCParticleClassifier.MCPClass;
10
11 /**
12 * CreateFinalStateMCParticleList writes a list of Final state MCParticles to event. The choices are GeneratorFS(Gen) and SimulatorFS(Sim). The final state particles are those at the end of the
13 * decay/interaction chain for the choice(Gen or Sim). By default no particles created without destroying the parent are included. This can be overridden by the setKeepContinuousXXX methods. Decays
14 * and interactions beyond a radius or z position can be ignored by using the setRadiusCut and setZCut methods. The collection name of the final state particles can also be set.
15 * @author Ron Cassell
16 */
17 public class CreateFinalStateMCParticleList extends Driver {
18 private MCPClass fs;
19 private String FStype;
20 private String CollectionName;
21 private double Rcut;
22 private double Zcut;
23 private boolean keepce;
24 private boolean keepcp;
25 private boolean keepch;
26 private final static String[] types = { "Gen", "Sim" };
27 private int itype;
28 private MCParticleClassifier cl;
29
30 /**
31 * Set up the defaults for creating the lists. Type parameter must be Gen or Sim.
32 *
33 * @param type
34 */
35 public CreateFinalStateMCParticleList(String type) {
36 fs = MCPClass.GEN_FINAL_STATE;
37 cl = new MCParticleClassifier();
38 FStype = type;
39 if (FStype.compareTo(types[0]) == 0) {
40 itype = 0;
41 } else if (FStype.compareTo(types[1]) == 0) {
42 itype = 1;
43 } else {
44 itype = 0;
45 System.out.println("CreateFinalStateMCParticleList created with invalid type " + FStype + ": Defaulting to " + types[0]);
46 FStype = types[0];
47 }
48 Rcut = 99999.;
49 Zcut = 99999.;
50 keepce = false;
51 keepcp = false;
52 keepch = false;
53 CollectionName = FStype + "FinalStateParticles";
54 }
55
56 /**
57 * Set the output collection name
58 *
59 * @param name
60 */
61 public void setCollectionName(String name) {
62 CollectionName = name;
63 }
64
65 /**
66 * Set a radius cut beyond which interactions and decays are ignored.
67 *
68 * @param rc
69 */
70 public void setRadiusCut(double rc) {
71 Rcut = rc;
72 }
73
74 /**
75 * Set a z cut beyond which interactions and decays are ignored.
76 *
77 * @param zc
78 */
79 public void setZCut(double zc) {
80 Zcut = zc;
81 }
82
83 /**
84 * Keep electrons created without destroying parent. (deltas, comptons)
85 */
86 public void setKeepContinuousElectrons() {
87 keepce = true;
88 }
89
90 /**
91 * Keep photons created without destroying parent. (brem, nuclear elastic and pseudo-elastic)
92 */
93 public void setKeepContinuousPhotons() {
94 keepcp = true;
95 }
96
97 /**
98 * Keep hadrons created without destroying parent. (nuclear elastic and pseudo-elastic)
99 */
100 public void setKeepContinuousHadrons() {
101 keepch = true;
102 }
103
104 /**
105 * Process and event. Create the final state particle list and store it in event.
106 *
107 * @param event
108 */
109 public void process(EventHeader event) {
110 //
111 // Create a list to hold the final state particles
112 //
113 List<MCParticle> fslist = new ArrayList<MCParticle>();
114 //
115 // Get the list of all MCParticles
116 //
117 List<MCParticle> all = (List<MCParticle>) event.get("MCParticle");
118 //
119 // If Generator final state particles requested, use the MCParticleClassifier
120 // to identify the final state particles, and add them to the list
121 //
122 if (itype == 0) {
123 boolean simulated = false;
124 for (MCParticle p : all) {
125 if (p.getSimulatorStatus().isDecayedInTracker() || p.getSimulatorStatus().isDecayedInCalorimeter() || p.getSimulatorStatus().hasLeftDetector() || p.getSimulatorStatus().isStopped()) {
126 simulated = true;
127 continue;
128 }
129 }
130 if (!simulated) {
131 for (MCParticle p : all) {
132 if (p.getGeneratorStatus() == MCParticle.FINAL_STATE)
133 fslist.add(p);
134 }
135 } else {
136 for (MCParticle p : all) {
137 if (cl.getClassification(p) == fs) {
138 fslist.add(p);
139 }
140 }
141 }
142 }
143 //
144 // Simulator final state particles requested
145 //
146 else {
147 //
148 // Check if we are keeping any continuous process created particles
149 //
150 boolean keepsome = keepce || keepcp || keepch;
151 //
152 // Loop over all MCParticles
153 //
154 for (MCParticle p : all) {
155 //
156 // Never keep backscatter particles
157 //
158 if (p.getSimulatorStatus().isBackscatter())
159 continue;
160 //
161 // Don't keep particles the Simulator never saw
162 //
163 boolean inSim = p.getSimulatorStatus().isDecayedInTracker() || p.getSimulatorStatus().isDecayedInCalorimeter() || p.getSimulatorStatus().hasLeftDetector() || p.getSimulatorStatus().isStopped();
164 if (!inSim)
165 continue;
166 //
167 // Find out if the particle has endpoint daughters
168 //
169 boolean hasepd = false;
170 for (MCParticle d : p.getDaughters()) {
171 if (d.getSimulatorStatus().isBackscatter())
172 continue;
173 if (d.getSimulatorStatus().vertexIsNotEndpointOfParent())
174 continue;
175 hasepd = true;
176 break;
177 }
178 //
179 // If the particle has endpoint daughters it is not final state
180 //
181 if (hasepd)
182 continue;
183 //
184 // This is a simulator final state particle. If it is a result of a
185 // continuous process, check to see if we keep it.
186 //
187 MCParticle pp = getFirstContinuousParticle(p);
188 if (pp != null) {
189 if (keepsome) {
190 if (keepThis(p, pp))
191 fslist.add(p);
192 }
193 }
194 //
195 // Add the final state particle to the list
196 //
197 else {
198 fslist.add(p);
199 }
200 }
201 }
202 //
203 // We now have the complete list of final state particles. Check for
204 // production vertex cuts
205 //
206 if ((Rcut < 9999.) || (Zcut < 9999.)) {
207 //
208 // Create a remove list and an add list
209 //
210 List<MCParticle> removelist = new ArrayList<MCParticle>();
211 List<MCParticle> addtolist = new ArrayList<MCParticle>();
212 //
213 // Loop over the final state list
214 //
215 for (MCParticle p : fslist) {
216 //
217 // Check vertex against cuts
218 //
219 Hep3Vector vtx = p.getOrigin();
220 if ((Math.abs(vtx.z()) > Zcut) || (Math.sqrt(vtx.x() * vtx.x() + vtx.y() * vtx.y()) > Rcut)) {
221 //
222 // Vertex failed cuts. Remove it.
223 //
224 removelist.add(p);
225 //
226 // Unless a continuous process particle, add the parent to the list of particles to add back in
227 //
228 if (!p.getSimulatorStatus().vertexIsNotEndpointOfParent()) {
229 if (!addtolist.contains(p.getParents().get(0)))
230 addtolist.add(p.getParents().get(0));
231 }
232 }
233 }
234 //
235 // Remove all the particles in the remove list from the final state list
236 //
237 for (MCParticle p : removelist) {
238 fslist.remove(p);
239 }
240 //
241 // Loop over the add list
242 //
243 for (MCParticle p : addtolist) {
244 Hep3Vector vtx = p.getOrigin();
245 MCParticle pp = p;
246 boolean hasparent = true;
247 //
248 // Trace parentage until cut is passed (or run out of parents or break chain with continuous particle)
249 //
250 while ((hasparent) && ((Math.abs(vtx.z()) > Zcut) || (Math.sqrt(vtx.x() * vtx.x() + vtx.y() * vtx.y()) > Rcut))) {
251 if (pp.getSimulatorStatus().vertexIsNotEndpointOfParent()) {
252 pp = null;
253 } else {
254 pp = pp.getParents().get(0);
255 }
256 if (pp == null) {
257 hasparent = false;
258 } else {
259 vtx = pp.getOrigin();
260 }
261 }
262 //
263 // Add the particle to the fs list if it exists
264 //
265 if (!(pp == null)) {
266 if (!fslist.contains(pp))
267 fslist.add(pp);
268 }
269 }
270 }
271 //
272 // Write the collection to event
273 //
274 event.put(CollectionName, fslist);
275 event.getMetaData(fslist).setSubset(true);
276 }
277
278 /**
279 * Decide if a particle resulting from continuous production should be kept as a final state particle.
280 *
281 * @param p - particle in question
282 * @param pp - first continuous production particle in parentage chain of p
283 */
284 public boolean keepThis(MCParticle p, MCParticle pp) {
285 boolean keepit = false;
286 //
287 // if first vneop particle is the particle in question, only need to check it for validity
288 //
289 if (p == pp) {
290 if (Math.abs(p.getPDGID()) == 11) {
291 if (keepce)
292 keepit = true;
293 } else if (p.getPDGID() == 22) {
294 if (keepcp)
295 keepit = true;
296 } else {
297 if (keepch)
298 keepit = true;
299 }
300 return keepit;
301 }
302 //
303 // otherwise, check the whole chain for validity
304 //
305 else {
306 //
307 // First check the particle in question
308 //
309 if (p.getSimulatorStatus().vertexIsNotEndpointOfParent()) {
310 if (Math.abs(p.getPDGID()) == 11) {
311 if (keepce)
312 keepit = true;
313 } else if (p.getPDGID() == 22) {
314 if (keepcp)
315 keepit = true;
316 } else {
317 if (keepch)
318 keepit = true;
319 }
320 if (!keepit)
321 return keepit;
322 }
323 //
324 // Particle itself is valid, check the chain
325 //
326 keepit = true;
327 MCParticle ppp = p;
328 while ((keepit) && (ppp != pp)) {
329 ppp = ppp.getParents().get(0);
330 if (ppp.getSimulatorStatus().vertexIsNotEndpointOfParent()) {
331 keepit = false;
332 if (Math.abs(ppp.getPDGID()) == 11) {
333 if (keepce)
334 keepit = true;
335 } else if (ppp.getPDGID() == 22) {
336 if (keepcp)
337 keepit = true;
338 } else {
339 if (keepch)
340 keepit = true;
341 }
342 }
343 }
344 return keepit;
345 }
346 }
347
348 /**
349 * Find and return the top VNEOP particle in parentage chain
350 *
351 * @param p - particle in question
352 */
353 public MCParticle getFirstContinuousParticle(MCParticle p) {
354 MCParticle rp = null;
355 if (p.getSimulatorStatus().vertexIsNotEndpointOfParent())
356 rp = p;
357 MCParticle pp = p;
358 while (pp.getParents().size() == 1) {
359 pp = pp.getParents().get(0);
360 if (pp.getSimulatorStatus().vertexIsNotEndpointOfParent())
361 rp = pp;
362 }
363 return rp;
364 }
365 }