View Javadoc

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   * -Changed to allow no detector.  By default, now uses a B field of (0,0,0).
25   * -Removed unused trackingZMax variable.  FIXME: Should be using zmax?
26   * -Removed unused flags variable.  FIXME: Should be using collection flags?
27   *
28   * @author tonyj
29   * @version $Id: MCParticleConverter.java,v 1.13 2011/11/14 23:02:14 jeremy Exp $
30   */
31  class MCParticleConverter implements HepRepCollectionConverter
32  {
33  	private static final double[] IP = { 0,0,0 };
34  	//private boolean _noDetector = false;
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  		//int flags = meta.getFlags();
45  		Detector detector = null;
46  		/*  try
47          {
48  		 */
49  		try {
50  			detector = event.getDetector();
51  		}
52  		catch (Exception x)
53  		{}
54  		
55  		double trackingRMax = 10000;
56  		//double trackingZMax = 20000;
57  
58  		if (detector != null) {
59  			trackingRMax = detector.getConstants().get("tracking_region_radius").getValue();
60  			//trackingZMax = detector.getConstants().get("tracking_region_zmax").getValue();
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 					// Workaround for simdet
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 					// Use the helix swimmer to swim to end of tracking region
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 						//TODO are there nuetral types other than photon, neutrino, and neutral hadron?
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 //					HepRepInstance instanceX = factory.createHepRepInstance(charge == 0 ? neutral : charged, charge == 0 ? neutralType : chargedType);
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 					// if particle starts at origin and has no apprecaible pT, don't draw
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; // 1cm
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 				// Just ignore it for now.
197 			}
198 		}
199 		/*
200         }
201         catch(Exception ex)
202         {
203             // Just ignore it for now.
204             if(!_noDetector)
205             {
206                 System.out.println(ex);
207                 System.out.println("Cannot display MCParticles without a detector!");
208                 _noDetector = true;
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 }