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.EventHeader;
9 import org.lcsim.geometry.IDDecoder;
10 import org.lcsim.recon.cat.util.Const;
11 import org.lcsim.util.Driver;
12 import org.lcsim.util.aida.AIDA;
13
14
15
16
17
18
19
20
21
22
23
24
25
26 final public class EmcalMipStubs extends Driver {
27
28 private AIDA aida = AIDA.defaultInstance();
29 private int debugLevel;
30 private int minEcalHits=2;
31 private int evtNo=0;
32 private boolean useSecondMipStub=false;
33 private boolean useThirdMipStub=false;
34 private double[] tempBase;
35 private double[] tempDir;
36 private double[] strtP;
37 private double[] mdlP;
38 private double[] endP;
39 private double lmax;
40 private double lmin;
41 private double deltaR;
42 private int minLayer=0;
43 private CircleFromPoints cir;
44 static final double k_cm= Const.cm;
45 static final double barrelRadius = 126.0 * k_cm;
46 static final double rotationalKappaCompensation = 0.0 / k_cm;
47 private int minLayerCut = 3;
48 boolean createStubTrackList = false;
49
50 static final double cosThetaEndcap = 0.8;
51 static final double alphaFactor = 0.015 / k_cm;
52 static String emcalMipStubsListName = "EmcalMipStubs";
53
54
55
56 public EmcalMipStubs(String listName, int debugL) {
57 emcalMipStubsListName = listName;
58 this.debugLevel = debugL;
59 if (debugLevel>=1) System.out.println("created EmcalMipStubs object");
60 this.tempBase = new double[]{0.,0.,0.};
61 this.tempDir = new double[]{0.,0.,0.};
62 this.strtP = new double[]{0.,0.,0.};
63 this.mdlP = new double[]{0.,0.,0.};
64 this.endP = new double[]{0.,0.,0.};
65 this.cir = new CircleFromPoints();
66 }
67
68 public void setMinEcalHits(int i){ minEcalHits =i;}
69
70 public void setDebugLevel(int i){ debugLevel = i;}
71
72
73
74
75 public boolean isMipStub(EventHeader event, Cluster c){
76 if (minLayer > minLayerCut) return false;
77 if (c.getCalorimeterHits().size()<minEcalHits) return false;
78 return true;
79 }
80
81 public boolean isEndcap(Cluster c){
82
83
84 double[] p=c.getPosition();
85 return (Math.abs(p[2])/Math.sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2])>0.8);
86 }
87
88
89
90
91
92 public void process( EventHeader event) {
93 evtNo++;
94
95 List distClusters = new ArrayList();
96 List<MipStub> listStubs = new ArrayList();
97 int cnum=0;
98 List<List<Cluster>> clusterSets = event.get(Cluster.class);
99 List<CalorimeterHit> hitsBar = event.get(CalorimeterHit.class,"EcalBarrHits");
100 IDDecoder theCell = event.getMetaData(hitsBar).getIDDecoder();
101
102 for (List<Cluster> clusters : clusterSets) {
103
104 String name = event.getMetaData(clusters).getName();
105 if (name.indexOf("Ecal")<0) continue;
106 if (debugLevel >=1) System.out.println("name of Calorimeter"+name);
107
108
109 for (Cluster c : clusters){
110 cnum++;
111 List cHits = c.getCalorimeterHits();
112 MipStub stub1=null,stub2=null,stub3=null;
113 setStartEndPoint(theCell, c, isEndcap(c));
114 if (debugLevel>=1){
115 }
116 if (isMipStub(event,c)){
117 stub1 = createMipStubFromFit(theCell, c);
118
119 if (useSecondMipStub) stub2 = createMipStubStraight(theCell, c);
120 else stub2=null;
121
122 if (useThirdMipStub || isEndcap(c)) stub3 = createMipStub(theCell, c);
123 else stub3=null;
124
125 if (stub1 != null){
126
127 listStubs.add(stub1);
128 if (debugLevel>=4 && stub1 != null ) stub1.debug();
129 }
130 if (stub2 != null){
131
132 listStubs.add(stub2);
133 if (debugLevel>=4 && stub2 != null ) stub2.debug();
134 }
135 if (stub3 != null){
136
137 listStubs.add(stub3);
138 if (debugLevel>=4 && stub3 != null ) stub3.debug();
139 }
140 }
141 }
142 }
143 if (debugLevel >=3) System.out.println("EmcalMipStubs n="+listStubs.size());
144 event.put(emcalMipStubsListName,listStubs);
145
146 if (createStubTrackList){
147 List<GarfieldTrack> listStubTracks = new ArrayList<GarfieldTrack>();
148 for (MipStub mStub : listStubs) {
149 GarfieldTrack newTrack = new GarfieldTrack(mStub, debugLevel);
150 newTrack.setTrackParameters();
151 listStubTracks.add(newTrack);
152 }
153 event.put("GarfieldStubTracks",listStubTracks);
154 }
155 }
156
157
158 public MipStub createUIowaMipstub(IDDecoder theCell, Cluster c) {
159
160 return createMipStubSimple(theCell, c);
161 }
162
163 public MipStub createMipStubFromFit(IDDecoder theCell, Cluster c) {
164 int i;
165 int nHit = c.getCalorimeterHits().size();
166 cir.calculate(0.,0.,strtP[0],strtP[1],endP[0],endP[1]);
167 GarfieldTrack newTrk = new GarfieldTrack();
168 newTrk.setDebugLevel(debugLevel-2);
169 double err = 0.2 * k_cm;
170 for (i=0;i<=2;i++) tempBase[i]=0;
171 newTrk.addHit( new GarfieldHit(endP,err,4,4));
172 newTrk.addHit( new GarfieldHit(mdlP,5.*err,3,3));
173 newTrk.addHit( new GarfieldHit(strtP,err,2,2));
174 newTrk.addHit( new GarfieldHit(tempBase,10.*err,1,1));
175
176 newTrk.setHelix(strtP,endP);
177 newTrk.hel.setKappa(cir.getKappa());
178 newTrk.calculateChi2();
179
180 newTrk.fullChi2Fit(0.1,2);
181 if (strtP[0]*strtP[0]+strtP[1]*strtP[1]>2500.*k_cm*k_cm){
182 newTrk.hel.setPointOnHelixWithXY(strtP[0],strtP[1]);} else {newTrk.hel.setPointOnHelixWithZ(strtP[2]);}
183
184 for (i=0;i<=2;i++){
185 tempBase[i]=newTrk.hel.getPointOnHelix(i);
186 tempDir[i]=newTrk.hel.dirAtPoint(i);
187 }
188 return new MipStub(tempBase,tempDir,newTrk.hel.kappa(),nHit, isEndcap(c),this.minLayer, c);
189 }
190
191 public MipStub createMipStubKyoko(IDDecoder theCell, Cluster c) {
192 double[] base = new double[]{0,0,0};
193 double[] dir = new double[]{0,0,0};
194 double[] ddir = new double[]{0,0,0};
195 double kappa = 0.;
196 int nHit = c.getCalorimeterHits().size();
197 int iCount = 0;
198 double[] p;
199 double Rmax = -10.;
200 double RRmax = -10.;
201 double Rmin = 1000000.;
202 double R;
203 List<CalorimeterHit> cHits = c.getCalorimeterHits();
204 int n = 0;
205 double sigmaxy = 0.;
206 double sigmax = 0.;
207 double sigmay =0.;
208 double sigmaxx =0.;
209 double sigmayy = 0.;
210 double sigmaxx2;
211 double[] b = new double[100];
212 double[] rr = new double[100];
213 double [] x = new double[100];
214 double [] y = new double[100];
215 double [] z = new double[100];
216 double avex, avey;
217 double ssxx, ssyy,ssxy;
218
219 for (CalorimeterHit ch: cHits) {
220 iCount++;
221 theCell.setID(ch.getCellID());
222 p = ch.getPosition();
223 R = Math.sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
224 if (debugLevel>=5) System.out.println("EmcalMipStub pos="+p[0]+" "+p[1]+" "+p[2]);
225 if (R>RRmax){
226 RRmax = R;
227 ddir=p;
228 }
229 if (R<Rmin){
230 Rmin = R;
231 base=p;
232 }
233
234
235
236
237
238 n = iCount;
239 x[n] = p[0];
240 y[n] = p[1];
241 z[n] = p[2];
242 }
243 dir[0]=ddir[0]-base[0];
244 dir[1]=ddir[1]-base[1];
245 dir[2]=ddir[2]-base[2];
246 double cotTheta= dir[2]/Math.sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
247
248 sort(x,y,z,n);
249
250 for(int i=1;i<=n;i++){
251 sigmaxy += x[i]*y[i];
252 sigmax += x[i];
253 sigmay += y[i];
254 sigmaxx += x[i]*x[i];
255 sigmayy += y[i]*y[i];
256 sigmaxx2 = sigmax*sigmax;
257
258 avex = sigmax/i;
259 avey = sigmay/i;
260 ssxx = sigmaxx-i*avex*avex;
261 ssyy = sigmayy-i*avey*avey;
262 ssxy = sigmaxy-i*avex*avey;
263 if(ssxx==0){
264 b[i] = 1.E15;
265 rr[i] = 1.;
266 } else if(ssyy==0){
267 b[i] = (i*sigmaxy-sigmax*sigmay)/(i*sigmaxx-sigmaxx2);
268 rr[i] = 1.;
269 } else{
270 b[i] = (i*sigmaxy-sigmax*sigmay)/(i*sigmaxx-sigmaxx2);
271 rr[i] = ssxy*ssxy/ssxx/ssyy;
272 }
273
274 R = Math.sqrt(x[i]*x[i]+y[i]*y[i]+z[i]*z[i]);
275 if (i<6 && R>Rmax){
276 Rmax = R;
277 dir[0] = x[i];
278 dir[1] = y[i];
279 dir[2] = z[i];
280 } else{
281 if(rr[i-1]>rr[i]){
282 rr[i] = rr[i-1];
283 b[i] = b[i-1];
284 x[i] = x[i-1];
285 y[i] = y[i-1];
286 z[i] = z[i-1];
287 }
288 dir[0] = x[i];
289 dir[1] = y[i];
290 dir[2] = z[i];
291 }
292 }
293 System.out.println("The correlation coefficient is "+rr[n]);
294
295 R = norm(dir);
296 if(R>Rmax) Rmax = R;
297
298
299
300 if(rr[n]>0.7){
301
302 dir[0] = dir[0] - base[0];
303 dir[1] = dir[1] - base[1];
304 dir[2] = dir[2] - base[2];
305
306 if(dir[0]>0){
307 dir[0] = 1/Math.sqrt(1+b[n]*b[n]);
308 dir[1] = b[n]/Math.sqrt(1+b[n]*b[n]);
309 } else{
310 dir[0] = -1/Math.sqrt(1+b[n]*b[n]);
311 dir[1] = -b[n]/Math.sqrt(1+b[n]*b[n]);
312 }
313 dir[2] = cotTheta*Math.sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
314 } else{
315 dir[0] = ddir[0] - base[0];
316 dir[1] = ddir[1] - base[1];
317 dir[2] = ddir[2] - base[2];
318 }
319
320 double angle = getAngleBaseDir(base,dir);
321
322 if (!isEndcap(c) && Math.abs(angle)>Math.PI/2.){
323 kappa = 0.;
324 dir[0]=base[0];
325 dir[1]=base[1];
326 dir[2]=base[2];
327 } else{
328 kappa= 0.02*angle;
329 double maxKappaValue =1.2 *(1./( 0.5 * Math.sqrt(base[0]*base[0]+base[1]*base[1])));
330 if (kappa < -maxKappaValue) kappa = -maxKappaValue;
331 if (kappa > maxKappaValue) kappa = maxKappaValue;
332 }
333 return new MipStub(base,dir,kappa,nHit, isEndcap(c), this.minLayer);
334 }
335
336
337
338 public MipStub createMipStub(IDDecoder theCell, Cluster c) {
339 double kappa = 0.;
340 int nHit = c.getCalorimeterHits().size();
341 int iCount = 0;
342 double[] p;
343 double lngth;
344 if (lmax-lmin<1.E-1){
345 return new MipStub(strtP,strtP,kappa,nHit, isEndcap(c), this.minLayer);
346 }
347
348 int nBase = 0;
349 int nEnd = 0;
350 tempBase[0]=0.;
351 tempBase[1]=0.;
352 tempBase[2]=0.;
353 tempDir[0]=0.;
354 tempDir[1]=0.;
355 tempDir[2]=0.;
356 List<CalorimeterHit> cHits = c.getCalorimeterHits();
357 for (CalorimeterHit ch: cHits) {
358 theCell.setID(ch.getCellID());
359 p = ch.getPosition();
360 if (isEndcap(c)) lngth = Math.abs(p[2]);
361 else lngth = Math.sqrt(p[0]*p[0]+p[1]*p[1]);
362
363 if (lngth<lmin+ 2. * k_cm) {
364 nBase++;
365 tempBase[0]+=p[0];
366 tempBase[1]+=p[1];
367 tempBase[2]+=p[2];
368 } else if (lngth < lmin + 0.5 * (lmax-lmin)) {
369 nEnd++;
370 tempDir[0]+=p[0];
371 tempDir[1]+=p[1];
372 tempDir[2]+=p[2];
373 }
374 }
375
376 if (nEnd==0 || nBase==0){
377 tempBase[0]=strtP[0];
378 tempBase[1]=strtP[1];
379 tempBase[2]=strtP[2];
380 tempDir[0]=endP[0]-strtP[0];
381 tempDir[1]=endP[1]-strtP[1];
382 tempDir[2]=endP[2]-strtP[2];
383 } else{
384 tempBase[0]=tempBase[0]/(double) nBase;
385 tempBase[1]=tempBase[1]/(double) nBase;
386 tempBase[2]=tempBase[2]/(double) nBase;
387 tempDir[0]=(tempDir[0]/(double) nEnd)-tempBase[0];
388 tempDir[1]=(tempDir[1]/(double) nEnd)-tempBase[1];
389 tempDir[2]=(tempDir[2]/(double) nEnd)-tempBase[2];
390 }
391 double angle = getAngleBaseDir(tempBase,tempDir);
392 if (!isEndcap(c) && Math.abs(angle)>Math.PI/2.){
393 kappa = 0.;
394 tempDir[0]=tempBase[0];
395 tempDir[1]=tempBase[1];
396 tempDir[2]=tempBase[2];
397 } else{
398 kappa= alphaFactor*angle*barrelRadius/Math.sqrt(tempBase[0]*tempBase[0]+
399 tempBase[1]*tempBase[1]);
400
401 double maxKappaValue =1.2 *(1./( 0.5 * Math.sqrt(tempBase[0]*tempBase[0]+
402 tempBase[1]*tempBase[1])));
403 if (kappa < -maxKappaValue) kappa = -maxKappaValue;
404 if (kappa > maxKappaValue) kappa = maxKappaValue;
405 double phiRot = -rotationalKappaCompensation*kappa * Math.abs(deltaR);
406 double dirx=Math.cos(phiRot)*tempDir[0]-Math.sin(phiRot)*tempDir[1];
407 double diry=Math.sin(phiRot)*tempDir[0]+Math.cos(phiRot)*tempDir[1];
408 tempDir[0]=dirx;
409 tempDir[1]=diry;
410 }
411 return new MipStub(tempBase,tempDir,kappa,nHit, isEndcap(c), this.minLayer);
412 }
413
414
415 public MipStub createMipStubSimple(IDDecoder theCell, Cluster c) {
416 double kappa = 0.;
417 int nHit = c.getCalorimeterHits().size();
418 int iCount = 0;
419 this.lmax = -10.;
420 this.lmin = 1000000.;
421 double lngth;
422 double rxy;
423 List<CalorimeterHit> cHits = c.getCalorimeterHits();
424 if (lmax-lmin<1.E-1){return new MipStub(strtP,strtP,0.,nHit,isEndcap(c), this.minLayer);} else{
425 tempDir[0]=endP[0]-strtP[0];
426 tempDir[1]=endP[1]-strtP[1];
427 tempDir[2]=endP[2]-strtP[2];
428 MipStub stub = new MipStub(strtP,tempDir,kappa,nHit,isEndcap(c), this.minLayer);
429 kappa = alphaFactor* stub.getAngleBaseDir()*barrelRadius/Math.sqrt(strtP[0]*strtP[0]+strtP[1]*strtP[1]);
430
431 double maxKappaValue =1.2 *(1./( 0.5 * Math.sqrt(strtP[0]*strtP[0]+
432 strtP[1]*strtP[1])));
433 if (kappa < -maxKappaValue) kappa = -maxKappaValue;
434 if (kappa > maxKappaValue) kappa = maxKappaValue;
435 return new MipStub(strtP,tempDir,kappa,nHit,isEndcap(c), this.minLayer);
436 }
437 }
438
439 public MipStub createMipStubStraight(IDDecoder theCell, Cluster c) {
440 double kappa = 0.;
441 int nHit = c.getCalorimeterHits().size();
442 int iCount = 0;
443 this.lmax = -10.;
444 this.lmin = 1000000.;
445 double lngth;
446 double rxy;
447 List<CalorimeterHit> cHits = c.getCalorimeterHits();
448 if (lmax-lmin<1.E-1){return null;} else{
449 tempDir[0]=endP[0]-strtP[0];
450 tempDir[1]=endP[1]-strtP[1];
451 tempDir[2]=endP[2]-strtP[2];
452 return new MipStub(strtP,tempDir,kappa,nHit,isEndcap(c), this.minLayer);
453 }
454 }
455
456 private void setStartEndPoint(IDDecoder theCell, Cluster c, boolean isEndcap){
457 int i;
458 List<CalorimeterHit> cHits = c.getCalorimeterHits();
459 double lngth;
460 double rxy;
461 double rxyMin=-10.;
462 double rxyMax=1000000;
463 this.minLayer = 1000000;
464 this.lmax= -1.E10;
465 this.lmin = 1.E10;
466 for (CalorimeterHit ch: cHits) {
467 int iCount = 0;
468 iCount++;
469 theCell.setID(ch.getCellID());
470 double[] p = ch.getPosition();
471 if (isEndcap) lngth = Math.abs(p[2]);
472 else lngth = Math.sqrt(p[0]*p[0]+p[1]*p[1]);
473 rxy = Math.sqrt(p[0]*p[0]+p[1]*p[1]);
474 if (debugLevel>=5) System.out.println("EmcalMipStub pos="+p[0]+" "+p[1]+" "+p[2]);
475 if (lngth>lmax) {
476 this.lmax = lngth;
477 rxyMax = rxy;
478 for (i=0;i<=2;i++) endP[i]=p[i];
479 }
480 if (lngth<lmin) {
481 this.lmin = lngth;
482 rxyMin = rxy;
483 for (i=0;i<=2;i++) strtP[i]=p[i];
484 int theLayer = theCell.getLayer();
485 if (theLayer < minLayer) this.minLayer = theLayer;
486 }
487 }
488
489 this.deltaR = Math.abs(rxyMax - rxyMin);
490 for (i=0;i<=2;i++) mdlP[i]=0.;
491 int nmdl=0;
492 List<CalorimeterHit> cHits2 = c.getCalorimeterHits();
493 for (CalorimeterHit ch: cHits2) {
494
495
496 theCell.setID(ch.getCellID());
497 double[] p = ch.getPosition();
498 if (isEndcap) lngth = Math.abs(p[2]);
499 else lngth = Math.sqrt(p[0]*p[0]+p[1]*p[1]);
500 if (Math.abs(lngth-0.5*(lmax+lmin))<0.15*(lmax-lmin)){
501 for (i=0;i<=2;i++) mdlP[i]=mdlP[i]+p[i];
502 nmdl++;
503 }
504 }
505 if (norm(mdlP)<1.*k_cm){
506 for (i=0;i<=2;i++) mdlP[i]=0.5*(strtP[i]+endP[i]);} else{ for (i=0;i<=2;i++) mdlP[i]=mdlP[i]/nmdl;}
507 }
508
509
510 private double norm(double[]x){
511 return Math.sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]);
512 }
513
514 public static void sort(double[]x,double[]y,double[]z, int n){
515 double[] RR = new double[100];
516 for(int i=1;i<=n;i++){
517 RR[i] = Math.sqrt(x[i]*x[i]+y[i]*y[i]+z[i]*z[i]);
518 }
519 for(int i=1;i<=n;i++){
520 int min = i;
521 for(int j=i;j<=n;j++){
522 if(RR[j]<RR[i]) min=j;
523 }
524 double tmpx,tmpy,tmpz;
525 tmpx = x[i];
526 tmpy = y[i];
527 tmpz = z[i];
528 x[i] = x[min];
529 y[i] = y[min];
530 z[i] = z[min];
531 x[min] = tmpx;
532 y[min] = tmpy;
533 z[min] = tmpz;
534 }
535 }
536
537 private double getAngleBaseDir(double[] base, double[] dir) {
538 double angleBaseDir;
539 double lBase = Math.sqrt(base[0]*base[0]+base[1]*base[1]);
540 double lDir = Math.sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
541 if (lBase*lDir>1.E-9){
542 double c=(base[0]*dir[0]+base[1]*dir[1])/(lBase*lDir);
543 if (Math.abs(c)>1.){ angleBaseDir = 0.;} else {angleBaseDir = Math.acos(c);}
544 if (base[0]*dir[1]-base[1]*dir[0]<0.) angleBaseDir=-angleBaseDir;
545 } else angleBaseDir = 0.;
546 return angleBaseDir;
547 }
548
549 public void setMinLayerCut(int val){minLayerCut =val;}
550 public void setCreateStubTrackList(boolean val){createStubTrackList =val;}
551 public void setUseSecondMipStub(boolean val){useSecondMipStub = val;}
552 public void setUseThirdMipStub(boolean val){useThirdMipStub = val;}
553
554 }
555
556
557