1 package org.lcsim.recon.tracking.trfcyl;
2
3 import org.lcsim.recon.tracking.trfbase.Interactor;
4 import org.lcsim.recon.tracking.trfbase.ETrack;
5 import org.lcsim.recon.tracking.trfbase.TrackVector;
6 import org.lcsim.recon.tracking.trfbase.TrackError;
7 import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
8
9 /**
10 * This class modifies the covariance matrix of a track which
11 * has been propagated to an InteractingLayer containing
12 * a LayerCylinder. The modifications correspond to multiple
13 *scattering in a thick cylindrical shell whose material is
14 *represented by the number of radiation lengths.
15 *
16 *
17 *@author Norman A. Graf
18 *@version 1.0
19 *
20 */
21
22 public class ThickCylMs extends Interactor
23 {
24
25 private double _pathLength; // Actual path length in cm
26 private double _radLength; // Radiation length of this cylinder
27
28 /**
29 *
30 * Construct an instance from the number of radiation
31 * lengths of the cylindrical shell material.
32 * The Interactor is constructed with the
33 * appropriate number of radiation lengths.
34 *
35 * @param pathLength The physical thickness of this cylindrical shell (in cm).
36 * @param radLength The thickness of the material in radiation lengths.
37 */
38
39 public ThickCylMs(double pathLength, double radLength)
40 {
41 _pathLength = pathLength;
42 _radLength = radLength;
43 }
44
45
46 /**
47 *Make a clone of this object.
48 *
49 * @return A Clone of this instance.
50 */
51 public Interactor newCopy()
52 {
53 return new ThickCylMs(_pathLength, _radLength);
54 }
55
56
57 /**
58 *Interact the given track in this cylindrical shell,
59 *using the thick material approximation for multiple scattering.
60 *Note that both the track parameters and the
61 *covariance matrix are updated to reflect the uncertainty caused
62 *by traversing the cylindrical shell of material.
63 *
64 * @param tre The ETrack to scatter.
65 */
66 public void interact(ETrack tre)
67 {
68
69 //if called with pl = 0, get nasty error matrix (log 0)!
70 if(_pathLength==0 || _radLength==0) return;
71
72
73 // This can only be used with cylinders... check that we have one..
74 // SurfCylinder cyl(10.0); //cylinder with radius 10 cm
75 // SurfacePtr srf = tre.get_surface();
76 // assert( srf->get_pure_type() == cyl.get_pure_type() );
77
78
79 TrackError cleanError = tre.error();
80 TrackError newError = cleanError;
81
82 TrackVector theVec = tre.vector();
83
84 //calculate q/p from q/pt
85 double trackMomentum = Math.abs(theVec.get(SurfCylinder.IQPT)*Math.cos(Math.atan(theVec.get(SurfCylinder.ITLM))));
86
87
88 //calculate dE/dx = a * p * s / X (a from fit to mc):
89 double dE = 0.01*(_pathLength/_radLength)/(trackMomentum);
90
91
92 //set the rms scattering appropriate to this momentum
93 //Theta = (0.0136 GeV)*(z/p)*(sqrt(thickness/rad_length))*(1+0.038*log(thickness/rad_length))
94
95 //calculate rms scattering angle
96 double stdTheta = (0.0136)*trackMomentum*Math.sqrt(_pathLength/_radLength)*
97 (1 + 0.038*Math.log(_pathLength/_radLength));
98
99 //the rms scattering displacement, y^2:
100 // double stdY = (1/(sqrt(3.0)))*_pathLength*stdTheta;
101
102 //there is a correlation between y and theta
103 double correl = Math.sqrt(3.0)/2;
104
105
106 // The MS covariance matrix can now be set.
107 // **************** code for matrix ***********************//
108
109 double ThetaSqr = stdTheta*stdTheta;
110 double tlamSqr = theVec.get(SurfCylinder.ITLM)*theVec.get(SurfCylinder.ITLM);
111 //get the final radius from the track
112
113 double radius = tre.spacePoint().rxy();
114
115 double sSqr = _pathLength*_pathLength;
116 double rSqr = radius * radius;
117
118 double val = newError.get(SurfCylinder.IZ,SurfCylinder.IZ) + sSqr*ThetaSqr*(1 + tlamSqr)/3;
119 newError.set(SurfCylinder.IZ,SurfCylinder.IZ, val);
120 val = newError.get(SurfCylinder.IALF,SurfCylinder.IALF) + ThetaSqr*( (1 + tlamSqr) + sSqr/(3*rSqr) );;
121 newError.set(SurfCylinder.IALF,SurfCylinder.IALF,val);
122 // - 2*correl*sqrt((1 + tlamSqr)*sSqr/(3*rSqr)) );
123 val = newError.get(SurfCylinder.IPHI,SurfCylinder.IPHI) + sSqr*ThetaSqr/(3*rSqr);
124 newError.set(SurfCylinder.IPHI,SurfCylinder.IPHI,val);
125
126 //same as for thin cylinder MS
127 val = newError.get(SurfCylinder.ITLM,SurfCylinder.ITLM) + ThetaSqr*(1 + tlamSqr)*(1 + tlamSqr);
128 newError.set(SurfCylinder.ITLM,SurfCylinder.ITLM,val);
129 val = newError.get(SurfCylinder.IQPT,SurfCylinder.IQPT) + ThetaSqr*theVec.get(SurfCylinder.IQPT)*theVec.get(SurfCylinder.IQPT)*tlamSqr;
130 newError.set(SurfCylinder.IQPT,SurfCylinder.IQPT,val);
131
132
133 //Sort out the correlations
134 //Insert values for off diagonals & use symmetry to set lower.
135 val = newError.get(SurfCylinder.IZ,SurfCylinder.ITLM) + correl*Math.sqrt(ThetaSqr)*(1 + tlamSqr) * Math.sqrt( sSqr*ThetaSqr*(1 + tlamSqr)/3 );
136 newError.set(SurfCylinder.IZ,SurfCylinder.ITLM,val);
137 newError.set(SurfCylinder.ITLM,SurfCylinder.IZ,val);
138
139 val = newError.get(SurfCylinder.ITLM,SurfCylinder.IQPT) + ThetaSqr*theVec.get(SurfCylinder.ITLM)*theVec.get(SurfCylinder.IQPT)*(1.0+tlamSqr);
140 newError.set(SurfCylinder.ITLM,SurfCylinder.IQPT,val);
141 newError.set(SurfCylinder.IQPT,SurfCylinder.ITLM, val);
142
143 val = newError.get(SurfCylinder.IALF,SurfCylinder.IPHI) + correl * ThetaSqr * Math.sqrt(sSqr/(3*rSqr) * ( (1 + tlamSqr) + sSqr/(3*rSqr)) );
144 newError.set(SurfCylinder.IALF,SurfCylinder.IPHI,val);
145 newError.set(SurfCylinder.IPHI,SurfCylinder.IALF,val);
146
147
148 // ********************************************************//
149
150
151 //correct for dE/dx:
152 // System.out.println("qoverpt: "+theVec.get(SurfCylinder.IQPT) +"; dE: "+dE);
153 // if ( theVec.get(SurfCylinder.IQPT)<0) theVec.set(SurfCylinder.IQPT, 1/( 1/(theVec.get(SurfCylinder.IQPT)) + dE ) );
154 // else theVec.set(SurfCylinder.IQPT, 1/( 1/(theVec.get(SurfCylinder.IQPT)) - dE ) );
155 // val = newError.get(SurfCylinder.IQPT,SurfCylinder.IQPT) + (theVec.get(SurfCylinder.IQPT)*theVec.get(SurfCylinder.IQPT))*0.5*( dE*dE );
156 // newError.set(SurfCylinder.IQPT,SurfCylinder.IQPT,val);
157
158 // tre.set_vector(theVec);
159
160 tre.setError( newError );
161 }
162
163 /**
164 *Return the number of radiation lengths.
165 *
166 * @return The thickness of the scattering material in radiation lengths.
167 */
168 public double radLength()
169 {
170 return _radLength;
171 }
172
173 /**
174 *Return the thickness of this cylindrical shell in cm.
175 *
176 * @return The thickness of the cylindrical shell in cm.
177 */
178 public double pathLength()
179 {
180 return _pathLength;
181 }
182
183
184 /**
185 *output stream
186 *
187 * @return A String representation of this instance.
188 */
189 public String toString()
190 {
191 return "ThickCylMs with "+_radLength+" radiation lengths and "+_pathLength+" thickness.";
192 }
193
194 }