1 /*
2 * CylEloss_Test.java
3 *
4 * Created on July 24, 2007, 8:55 PM
5 *
6 * $Id: CylEloss_Test.java,v 1.1.1.1 2010/04/08 20:38:00 jeremy Exp $
7 */
8
9 package org.lcsim.recon.tracking.trfcyl;
10
11 import junit.framework.TestCase;
12 import org.lcsim.recon.tracking.trfbase.ETrack;
13 import org.lcsim.recon.tracking.trfbase.Interactor;
14 import org.lcsim.recon.tracking.trfbase.PropDir;
15 import org.lcsim.recon.tracking.trfbase.Surface;
16 import org.lcsim.recon.tracking.trfbase.TrackError;
17 import org.lcsim.recon.tracking.trfbase.TrackVector;
18 import org.lcsim.recon.tracking.trfeloss.DeDx;
19 import org.lcsim.recon.tracking.trfeloss.DeDxFixed;
20 import org.lcsim.recon.tracking.trfutil.Assert;
21
22 import static java.lang.Math.cos;
23 import static java.lang.Math.abs;
24 import static java.lang.Math.atan;
25 import static java.lang.Math.sqrt;
26
27 /**
28 *
29 * @author Norman Graf
30 */
31 public class CylEloss_Test extends TestCase
32 {
33 private boolean debug;
34 /** Creates a new instance of CylEloss_Test */
35 public void testCylEloss()
36 {
37 String component = "CylEloss";
38 String ok_prefix = component + " (I): ";
39 String error_prefix = component + " test (E): ";
40
41 if(debug) System.out.println( ok_prefix
42 + "---------- Testing component " + component
43 + ". ----------" );
44
45 //********************************************************************
46 {
47 TrackVector trv = new TrackVector();
48 trv.set(0, 1.0);
49 trv.set(1, 1.0);
50 trv.set(2, 0.0);
51 trv.set(3, 0.0);
52 trv.set(4, 2.0);
53
54 double density = 1.0; // g/cm^3
55 double thickness = 2.; // cm
56
57 DeDx dedx = new DeDxFixed(density);
58 if(debug) System.out.println( ok_prefix + "Test constructor." );
59 CylEloss passIt = new CylEloss(thickness, dedx);
60 if(debug) System.out.println(passIt);
61 Assert.assertTrue(passIt.thickness() == thickness);
62 Assert.assertTrue(passIt.dEdX().equals(dedx));
63
64 TrackError initialError = new TrackError();
65
66 Surface srf = new SurfCylinder( 20.0 );
67
68 ETrack tmpTrack = new ETrack( srf, trv, initialError );
69
70 tmpTrack.setError( initialError );
71
72 TrackVector initialVector = tmpTrack.vector();
73 passIt.interact( tmpTrack );
74
75 TrackError finalError = tmpTrack.error();
76 TrackVector finalVector = tmpTrack.vector();
77
78 double particleMass = 0.13957;
79 double ptmax = 10000.; // GeV
80
81 double pinv = Math.abs(trv.get(4))*Math.cos(Math.atan(trv.get(3)));
82
83 // make sure pinv is greater than a threshold (1/ptmax)
84 // in this case assume q = 1, otherwise q = q/pt/abs(q/pt)
85
86 double sign = 1;
87 if(pinv < 1./ptmax)
88 pinv = 1./ptmax;
89 else
90 sign = trv.get(4)/Math.abs(trv.get(4));
91
92 double initialEnergy = Math.sqrt(1./pinv/pinv
93 +particleMass*particleMass);
94 double finalEnergy = initialEnergy;
95
96 double d = thickness/Math.abs(Math.cos(trv.get(2)))/Math.cos(Math.atan(trv.get(3)));
97
98 finalEnergy = dedx.loseEnergy(finalEnergy, d);
99 double sigmaEnergy = dedx.sigmaEnergy(initialEnergy, d);
100 double sigmaMomentum = sigmaEnergy;
101
102 if(finalEnergy<particleMass)finalEnergy=initialEnergy;
103
104 // now evaluate the final q/p and error(4,4)
105 double finalQoverP = sign/Math.sqrt(finalEnergy*finalEnergy-
106 particleMass*particleMass)/Math.cos(Math.atan(trv.get(3)));
107 double finalEr = sigmaMomentum*trv.get(4)*trv.get(4)*Math.cos(Math.atan(trv.get(3)));
108 finalEr *= finalEr;
109
110 if(debug) System.out.println("initial E= "+1/initialVector.get(4));
111 if(debug) System.out.println("interacted E= "+1/finalVector.get(4)+ " calculated E= "+1/finalQoverP);
112 Assert.assertTrue(Math.abs(finalVector.get(4)-finalQoverP)<0.00001);
113 Assert.assertTrue(Math.abs(finalError.get(4,4)-finalEr)<0.00001);
114
115 //********************************************************************
116 {
117 if(debug) System.out.println(ok_prefix + "Testing directionality");
118 ETrack tre1 = new ETrack(tmpTrack);
119 ETrack tre2 = new ETrack(tmpTrack);
120 if(debug) System.out.println("Track before interaction "+tmpTrack);
121 tre1.setForward();
122 tre1.setTrackForward();
123 passIt.interact(tre1);
124 if(debug) System.out.println("Track after forward interaction "+tre1);
125 //Should lose energy going forward
126 Assert.assertTrue(tre1.vector(4)>tmpTrack.vector(4));
127 tre2.setBackward();
128 tre2.setTrackBackward();
129 passIt.interact(tre2);
130 //Should gain energy going backward
131 Assert.assertTrue(tre2.vector(4)<tmpTrack.vector(4));
132 if(debug) System.out.println("Track after backward interaction "+tre2);
133 }
134 //********************************************************************
135 {
136 if(debug) System.out.println(ok_prefix + "Testing clone");
137 ETrack tre1 = new ETrack(tmpTrack);
138 ETrack tre2 = new ETrack(tmpTrack);
139 Assert.assertTrue(tre1.equals(tre2));
140 // clone
141 Interactor passIt2 = passIt.newCopy();
142 // assert they are not the same object
143 Assert.assertTrue( passIt2!=passIt );
144 // scatter both tracks
145 passIt2.interact(tre2);
146 passIt.interact(tre1);
147 // assert they're now equal
148 Assert.assertTrue( tre1.equals(tre2) );
149 // assert they're different from the original track
150 Assert.assertTrue( tre1.notEquals(tmpTrack) );
151 }
152 }
153 //********************************************************************
154
155 {
156 if(debug) System.out.println(ok_prefix + "Testing directionality");
157 TrackVector trv = new TrackVector();
158 trv.set(0, 1.0);
159 trv.set(1, 1.0);
160 trv.set(2, 0.0);
161 trv.set(3, 0.0);
162 trv.set(4, 0.0);
163
164 double density = 1.0; // g/cm^3
165 double thickness = 1.0; // cm
166
167 DeDx dedx = new DeDxFixed(density);
168 CylEloss passIt = new CylEloss(thickness, dedx);
169
170 TrackError initialError = new TrackError();
171
172 Surface srf = new SurfCylinder( 20.0 );
173
174 ETrack tmpTrack = new ETrack( srf, trv, initialError );
175
176 tmpTrack.setError( initialError );
177
178 TrackVector initialVector = tmpTrack.vector();
179 passIt.interact_dir( tmpTrack, PropDir.FORWARD );
180
181 TrackError finalError = tmpTrack.error();
182 TrackVector finalVector = tmpTrack.vector();
183
184 double particleMass = 0.13957;
185 double ptmax = 10000.; // GeV
186
187 double pinv = abs(trv.get(4))*cos(atan(trv.get(3)));
188
189 // make sure pinv is greater than a threshold (1/ptmax)
190 // in this case assume q = 1, otherwise q = q/pt/abs(q/pt)
191
192 int sign = 1;
193 if(pinv < 1./ptmax)
194 pinv = 1./ptmax;
195 else
196 sign = (int) (trv.get(4)/abs(trv.get(4)));
197
198 double initialEnergy = sqrt(1./pinv/pinv
199 +particleMass*particleMass);
200 double finalEnergy = initialEnergy;
201
202 double d = thickness/abs(cos(trv.get(2)))/cos(atan(trv.get(3)));
203
204 dedx.loseEnergy(finalEnergy, d);
205 double sigmaEnergy = dedx.sigmaEnergy(initialEnergy, d);
206 double sigmaMomentum = sigmaEnergy;
207
208 if(finalEnergy<particleMass)finalEnergy=initialEnergy;
209
210 // now evaluate the final q/p and error(4,4)
211 double finalQoverP = sign/sqrt(finalEnergy*finalEnergy-
212 particleMass*particleMass)/cos(atan(trv.get(3)));
213 double finalEr = sigmaMomentum*trv.get(4)*trv.get(4)*cos(atan(trv.get(3)));
214 finalEr *= finalEr;
215
216 Assert.assertTrue(abs(finalVector.get(4)-finalQoverP)<0.00001);
217 Assert.assertTrue(abs(finalError.get(4,4)-finalEr)<0.00001);
218
219 }
220
221
222 if(debug) System.out.println( ok_prefix
223 + "------------- All tests passed. -------------" );
224
225 }
226
227 }