1 /*
2 * XYPlaneEloss_Test.java
3 *
4 * Created on September 10, 2007, 11:33 AM
5 *
6 * $Id: XYPlaneEloss_Test.java,v 1.1.1.1 2010/04/08 20:38:00 jeremy Exp $
7 */
8
9 package org.lcsim.recon.tracking.trfxyp;
10
11 import junit.framework.TestCase;
12 import org.lcsim.recon.tracking.trfbase.ETrack;
13 import org.lcsim.recon.tracking.trfbase.PropDir;
14 import org.lcsim.recon.tracking.trfbase.Surface;
15 import org.lcsim.recon.tracking.trfbase.TrackError;
16 import org.lcsim.recon.tracking.trfbase.TrackVector;
17 import org.lcsim.recon.tracking.trfeloss.DeDx;
18 import org.lcsim.recon.tracking.trfeloss.DeDxFixed;
19 import org.lcsim.recon.tracking.trfzp.ZPlaneEloss;
20
21
22 import static java.lang.Math.abs;
23 import static java.lang.Math.sqrt;
24
25 /**
26 *
27 * @author Norman Graf
28 */
29 public class XYPlaneEloss_Test extends TestCase
30 {
31 private boolean debug;
32
33 /** Creates a new instance of XYPlaneEloss */
34 public void testXYPlaneEloss()
35 {
36 String component = "XYPlaneEloss";
37 String ok_prefix = component + " (I): ";
38 String error_prefix = component + " test (E): ";
39
40 if(debug) System.out.println( ok_prefix
41 + "---------- Testing component " + component
42 + ". ----------" );
43
44 //********************************************************************
45 TrackVector trv = new TrackVector();
46 trv.set(0,1.0);
47 trv.set(1, 1.0);
48 trv.set(2, 0.1);
49 trv.set(3, 0.2);
50 trv.set(4, 0.0);
51
52 double density = 1.0; // g/cm^3
53 double thickness = 1.0; // cm
54
55 DeDx dedx = new DeDxFixed(density);
56 XYPlaneEloss passIt = new XYPlaneEloss(thickness, dedx);
57
58 TrackError initialError = new TrackError();
59
60 Surface srf = new SurfXYPlane( 20.0, 1. );
61
62 ETrack tmpTrack = new ETrack( srf, trv, initialError );
63
64 tmpTrack.setError( initialError );
65
66 TrackVector initialVector = tmpTrack.vector();
67 passIt.interact_dir( tmpTrack, PropDir.FORWARD );
68
69 TrackError finalError = tmpTrack.error();
70 TrackVector finalVector = tmpTrack.vector();
71
72 double particleMass = 0.13957;
73 double ptmax = 10000.; // GeV
74
75 double pinv = abs(trv.get(4));
76
77 // make sure pinv is greater than a threshold (1/ptmax)
78 // in this case assume q = 1, otherwise q = q/pt/abs(q/pt)
79
80 int sign = 1;
81 if(pinv < 1./ptmax)
82 pinv = 1./ptmax;
83 else
84 sign = (int) (trv.get(4)/abs(trv.get(4)));
85
86 double initialEnergy = sqrt(1./pinv/pinv + particleMass*particleMass);
87 double finalEnergy = initialEnergy;
88
89 double d = thickness * sqrt(1. + trv.get(2)*trv.get(2) + trv.get(3)*trv.get(3));
90
91 dedx.loseEnergy(finalEnergy, d);
92 double sigmaEnergy = dedx.sigmaEnergy(initialEnergy, d);
93
94 if(finalEnergy<particleMass) finalEnergy=initialEnergy;
95
96 // now evaluate the final q/p and error(4,4)
97 double finalQoverP = sign/sqrt(finalEnergy*finalEnergy - particleMass*particleMass);
98 double finalEr = sigmaEnergy*trv.get(4)*trv.get(4);
99 finalEr *= finalEr;
100
101 assertTrue(abs(finalVector.get(4)-finalQoverP)<0.00001);
102 assertTrue(abs(finalError.get(4,4)-finalEr)<0.00001);
103
104 }
105 }