1 package org.lcsim.recon.cat;
2
3 import java.util.ArrayList;
4 import java.util.List;
5
6 import org.lcsim.event.CalorimeterHit;
7 import org.lcsim.event.Cluster;
8 import org.lcsim.event.MCParticle;
9
10
11
12
13
14
15
16
17
18
19
20
21 final public class MipStub implements Cluster {
22
23
24
25 public double[] base;
26 public double[] dir;
27 public double kappa;
28 private double baseError;
29 private double dirXYError;
30 private double dirZError;
31 private double kappaError;
32 private boolean isEndcap = false;
33 public double energy;
34 private double energyError;
35 public double angleBaseDir;
36
37 static final double k_cm = 10.;
38 public int debugLevel;
39 private int nHits=0;
40 private int minLayer=-1;
41 private MCParticle[] mcParticles = null;
42
43 private Cluster cluster = null;
44
45
46
47
48
49
50 public MipStub() {
51 this.base = new double[]{0,0,0};
52 this.dir = new double[]{0,0,0};
53 this.kappa = 0.;
54 this.energy = -1.;
55 this.angleBaseDir=0.;
56 nHits = 0;
57 normalizeDir();
58 setBaseError();
59 setDirError();
60 setKappaError();
61 }
62
63
64
65
66
67
68
69
70
71
72
73 public MipStub(double[] p, double[] d, double ka, int nHit, boolean isNdCap, int minLay) {
74 this.base = new double[]{p[0],p[1],p[2]};
75 this.dir = new double[]{d[0],d[1],d[2]};
76 this.kappa = ka;
77 this.energy = -1.;
78 this.nHits = nHit;
79 this.isEndcap = isNdCap;
80 this.minLayer = minLay;
81 normalizeDir();
82 setAngleBaseDir();
83 setBaseError();
84 setDirError();
85 setKappaError();
86 }
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101 public MipStub(double[] p, double[] d, double ka, int nHit, boolean isNdCap, int minLay, Cluster cluster) {
102 this.base = new double[]{p[0],p[1],p[2]};
103 this.dir = new double[]{d[0],d[1],d[2]};
104 this.kappa = ka;
105 this.energy = -1.;
106 this.nHits = nHit;
107 this.isEndcap = isNdCap;
108 this.minLayer = minLay;
109 normalizeDir();
110 setAngleBaseDir();
111 setBaseError();
112 setDirError();
113 setKappaError();
114 this.cluster = cluster;
115 }
116
117
118
119 public double kappa(){ return kappa;}
120 public int getnHits(){ return nHits;}
121 public double getPhi(){return Math.atan2(base[1],base[0]);}
122 public double getAngleBaseDir(){ return angleBaseDir;}
123 public double getDirXYError(){ return dirXYError;}
124 public double getDirZError(){ return dirZError;}
125 public double getBaseError(){return baseError;}
126 public double getKappaError(){return kappaError;}
127 public int getMinLayer(){return minLayer;}
128 public boolean isEndcap(){return isEndcap;}
129 public int charge() { return (kappa<0.) ? -1 : 1 ;}
130 public MCParticle[] getMCParticles() {return mcParticles;}
131
132 public double getEnergyError()
133 {
134 return energyError;
135 }
136
137 public void setEnergyError(double energyError)
138 {
139 this.energyError = energyError;
140 }
141
142
143
144 public void setMCParticles(MCParticle[] mcParticles) {this.mcParticles = mcParticles;}
145 public void setKappa(double ka){kappa=ka;}
146
147
148
149
150
151
152
153
154
155
156 private void normalizeDir(){
157
158 double d=Math.sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
159 if (d>1.E-20){
160 dir[0]=dir[0]/d;
161 dir[1]=dir[1]/d;
162 dir[2]=dir[2]/d;
163 } else {
164 dir[0]=1;
165 dir[1]=0.;
166 dir[2]=dir[2]*1.E10;
167 }
168 }
169
170
171
172
173 private void setAngleBaseDir() {
174 double lBase = Math.sqrt(base[0]*base[0]+base[1]*base[1]);
175 if (lBase>1.E-9){
176 double c=(base[0]*dir[0]+base[1]*dir[1])/(lBase);
177 if (Math.abs(c)>1.){ angleBaseDir = 0.;} else {angleBaseDir = Math.acos(c);}
178 if (base[0]*dir[1]-base[1]*dir[0]<0.) angleBaseDir=-angleBaseDir;
179 } else angleBaseDir = 0.;
180 }
181
182
183
184
185
186 private void setDirError(){
187 if (nHits>=10) dirXYError=20.;
188 else if (nHits >=5) dirXYError=25.;
189 else dirXYError= 45.;
190 dirXYError = Math.sqrt(dirXYError*dirXYError+900.*angleBaseDir*angleBaseDir);
191 dirZError = dirXYError / 5. * k_cm;
192
193 dirXYError = dirXYError*Math.PI/180.;
194 if (isEndcap()){
195 dirZError = dirZError * 10.;
196 dirXYError = dirXYError * 10.;
197 }
198 }
199
200 private void setBaseError(){ baseError = 1.0 * k_cm;}
201
202 private void setKappaError(){
203
204
205 if (nHits>=15) kappaError=0.03;
206 else if (nHits>=10) kappaError=0.035;
207 else if (nHits>=5) kappaError=0.05;
208 else kappaError=0.1;
209 if (isEndcap()) kappaError = kappaError * 10.;
210 }
211
212
213
214
215 void debug(){
216
217 System.out.println("MipStub base="+this.base[0]+" "+this.base[1]+" "+this.base[2]);
218 System.out.println("MipStub dir="+this.dir[0]+" "+this.dir[1]+" "+this.dir[2]+" kappa="+this.kappa+" angle="+angleBaseDir);
219 }
220
221
222
223
224
225
226 public List plot(){
227
228 List lst = new ArrayList();
229 double s = 2.*k_cm;
230 for (int j=0;j<=5;j++){
231 double[] pix=new double[]{0,0,0};
232 pix[0]=base[0]+s*dir[0]*j;
233 pix[1]=base[1]+s*dir[1]*j;
234 pix[2]=base[2]+s*dir[2]*j;
235 lst.add(pix);
236 }
237 GarfieldHelix hh = new GarfieldHelix(base,dir,kappa);
238 for (int j=0;j<=10;j++){
239 double[] pix=new double[]{0,0,0};
240 hh.setPointOnHelix(-8.*k_cm*j);
241 pix[0]=hh.getPointOnHelix(0);
242 pix[1]=hh.getPointOnHelix(1);
243 pix[2]=hh.getPointOnHelix(2);
244 lst.add(pix);
245 }
246 return lst;
247 }
248
249
250
251
252 public int getType() {return 0;}
253
254 public double[] getSubdetectorEnergies() {return new double[] {getEnergy()};}
255
256 public double[] getShape() {return new double[6];}
257
258
259 public double[] getPositionError() {return new double[] {baseError, baseError, baseError, baseError, baseError, baseError};}
260
261 public double[] getPosition() {return base;}
262
263 public double getITheta() {return 0.;}
264
265 public double getIPhi() {return 0.;}
266
267 public double[] getHitContributions() {
268 if (cluster == null) {
269 return new double[] {0.};
270 } else {
271 return cluster.getHitContributions();
272 }
273 }
274
275 public double getEnergy() {return energy;}
276
277 public double[] getDirectionError() {return new double[] {dirXYError, dirXYError, dirZError};}
278
279 public List<Cluster> getClusters() {
280 List<Cluster> list = new ArrayList<Cluster>(1);
281 if (cluster != null) {
282 list.add(cluster);
283 }
284 return list;
285 }
286
287 public List<CalorimeterHit> getCalorimeterHits() {
288 if (cluster == null) {
289 return new ArrayList<CalorimeterHit>(1);
290 } else {
291 return cluster.getCalorimeterHits();
292 }
293 }
294
295 public int getSize() {return nHits;}
296
297 public int getParticleId() {
298 return 0;
299 }
300
301 }
302