1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | #include <cassert> |
9 | #include <math.h> |
10 | #include <map> |
11 | |
12 | #include "BCAL/DBCALHit.h" |
13 | #include "BCAL/DBCALTDCHit.h" |
14 | #include "BCAL/DBCALPoint.h" |
15 | #include "BCAL/DBCALGeometry.h" |
16 | #include "BCAL/DBCALShower_factory_KLOE.h" |
17 | |
18 | #include "DANA/DApplication.h" |
19 | |
20 | #include "units.h" |
21 | |
22 | using namespace std; |
23 | |
24 | |
25 | |
26 | |
27 | jerror_t DBCALShower_factory_KLOE::init() |
28 | { |
29 | |
30 | ethr_cell=0.0001; |
31 | |
32 | CLUST_THRESH = 0.03; |
33 | |
34 | elyr = 1; |
35 | xlyr = 2; |
36 | ylyr = 3; |
37 | zlyr = 4; |
38 | tlyr = 5; |
39 | |
40 | |
41 | MERGE_THRESH_DIST = 40.0; |
42 | MERGE_THRESH_TIME = 2.5; |
43 | MERGE_THRESH_ZDIST = 30.0; |
44 | MERGE_THRESH_XYDIST = 40.0; |
45 | |
46 | |
47 | BREAK_THRESH_TRMS= 5.0; |
48 | |
49 | gPARMS->SetDefaultParameter( "BCALRECON:CLUST_THRESH", CLUST_THRESH ); |
50 | gPARMS->SetDefaultParameter( "BCALRECON:MERGE_THRESH_DIST", MERGE_THRESH_DIST ); |
51 | gPARMS->SetDefaultParameter( "BCALRECON:MERGE_THRESH_TIME", MERGE_THRESH_TIME ); |
52 | gPARMS->SetDefaultParameter( "BCALRECON:MERGE_THRESH_ZDIST", MERGE_THRESH_ZDIST ); |
53 | gPARMS->SetDefaultParameter( "BCALRECON:MERGE_THRESH_XYDIST", MERGE_THRESH_XYDIST ); |
54 | gPARMS->SetDefaultParameter( "BCALRECON:BREAK_THRESH_TRMS", BREAK_THRESH_TRMS ); |
55 | |
56 | |
57 | |
58 | |
59 | |
60 | |
61 | |
62 | |
63 | |
64 | |
65 | |
66 | |
67 | |
68 | |
69 | |
70 | |
71 | |
72 | |
73 | |
74 | |
75 | |
76 | |
77 | |
78 | m_scaleZ_p0 = 0.99798; |
79 | m_scaleZ_p1 = 0.000361096; |
80 | m_scaleZ_p2 = -2.17338e-06; |
81 | m_scaleZ_p3 = 1.32201e-09; |
82 | |
83 | m_nonlinZ_p0 = -0.0201272; |
84 | m_nonlinZ_p1 = 0.000103649; |
85 | m_nonlinZ_p2 = 0; |
86 | m_nonlinZ_p3 = 0; |
87 | |
88 | |
89 | return NOERROR; |
90 | } |
91 | |
92 | |
93 | |
94 | |
95 | jerror_t DBCALShower_factory_KLOE::brun(JEventLoop *loop, int32_t runnumber) |
96 | { |
97 | |
98 | DApplication* app = dynamic_cast<DApplication*>(loop->GetJApplication()); |
99 | DGeometry* geom = app->GetDGeometry(runnumber); |
100 | geom->GetTargetZ(m_z_target_center); |
101 | |
102 | vector<const DBCALGeometry*> bcalGeomVect; |
103 | loop->Get( bcalGeomVect ); |
104 | bcalGeom = bcalGeomVect[0]; |
105 | |
106 | |
107 | |
108 | |
109 | |
110 | ATTEN_LENGTH = bcalGeom->GetBCAL_attenutation_length(); |
111 | C_EFFECTIVE = bcalGeom->GetBCAL_c_effective(); |
112 | |
113 | fiberLength = bcalGeom->GetBCAL_length(); |
114 | zOffset = bcalGeom->GetBCAL_center(); |
115 | |
116 | |
117 | int modmin = 0; |
118 | int modmax = bcalGeom->GetBCAL_Nmodules(); |
119 | int rowmin1=0; |
120 | int rowmax1= bcalGeom->GetBCAL_NInnerLayers(); |
121 | int rowmin2= rowmax1; |
122 | int rowmax2= bcalGeom->GetBCAL_NOuterLayers()+rowmin2; |
123 | int colmin1=0; |
124 | int colmax1=bcalGeom->GetBCAL_NInnerSectors(); |
125 | int colmin2=0; |
126 | int colmax2=bcalGeom->GetBCAL_NOuterSectors(); |
127 | |
128 | float r_inner= bcalGeom->GetBCAL_inner_rad(); |
129 | |
130 | for (int i = (rowmin1+1); i < (rowmax1+1); i++){ |
131 | |
132 | int cellId = bcalGeom->cellId(1,i,1); |
133 | |
134 | rt[i]=bcalGeom->r(cellId)-r_inner; |
135 | } |
136 | |
137 | for (int i = (rowmin2+1); i < (rowmax2+1); i++){ |
138 | int cellId = bcalGeom->cellId(1,i,1); |
139 | rt[i]=bcalGeom->r(cellId)-r_inner; |
140 | } |
141 | |
142 | |
143 | float r[modulemax_bcal48][layermax_bcal10][colmax_bcal4]; |
144 | float phi[modulemax_bcal48][layermax_bcal10][colmax_bcal4]; |
145 | |
146 | |
147 | for (int k = modmin; k < modmax; k++){ |
148 | for (int i = rowmin1; i < rowmax1; i++){ |
149 | for (int j = colmin1; j < colmax1; j++){ |
150 | |
151 | |
152 | int cellId = bcalGeom->cellId(k+1,i+1,j+1); |
153 | r[k][i][j]=bcalGeom->r(cellId); |
154 | phi[k][i][j]=bcalGeom->phi(cellId); |
155 | |
156 | xx[k][i][j]=r[k][i][j]*cos(phi[k][i][j]); |
157 | yy[k][i][j]=r[k][i][j]*sin(phi[k][i][j]); |
158 | } |
159 | } |
160 | |
161 | for (int i = rowmin2; i < rowmax2; i++){ |
162 | for (int j = colmin2; j < colmax2; j++){ |
163 | int cellId = bcalGeom->cellId(k+1,i+1,j+1); |
164 | r[k][i][j]=bcalGeom->r(cellId); |
165 | phi[k][i][j]=bcalGeom->phi(cellId); |
166 | |
167 | xx[k][i][j]=r[k][i][j]*cos(phi[k][i][j]); |
168 | yy[k][i][j]=r[k][i][j]*sin(phi[k][i][j]); |
169 | } |
170 | } |
171 | } |
172 | |
173 | |
174 | |
175 | |
176 | |
177 | |
178 | return NOERROR; |
179 | } |
180 | |
181 | |
182 | |
183 | |
184 | jerror_t DBCALShower_factory_KLOE::evnt(JEventLoop *loop, uint64_t eventnumber) |
185 | { |
186 | |
187 | CellRecon(loop); |
188 | CeleToArray(); |
189 | PreCluster(loop); |
190 | ClusNorm(); |
191 | ClusAnalysis(); |
192 | Trakfit(); |
193 | |
194 | |
195 | vector<DBCALShower*> clusters; |
196 | int id = 0; |
197 | for (int i = 1; i < (clstot+1); i++){ |
198 | |
199 | int j=clspoi[i]; |
200 | |
201 | if( e_cls[j] < CLUST_THRESH ) continue; |
202 | |
203 | |
204 | DBCALShower *shower = new DBCALShower; |
205 | |
206 | |
207 | |
208 | |
209 | |
210 | |
211 | |
212 | |
213 | |
214 | |
215 | |
216 | |
217 | |
218 | |
219 | |
220 | |
221 | |
222 | |
223 | |
224 | |
225 | |
226 | |
227 | |
228 | |
229 | |
230 | |
231 | |
232 | |
233 | |
234 | vector<const DBCALPoint*> pointsInShower; |
235 | FindPointsInShower(j, loop, pointsInShower); |
236 | |
237 | |
238 | |
239 | |
240 | |
241 | |
242 | |
243 | |
244 | double E=0,x=0,y=0,z=0,t=0; |
245 | int N_cell=0; |
246 | double sig_x=0,sig_y=0,sig_z=0,sig_t=0; |
247 | double sum_z_wt=0; |
248 | for(unsigned int j=0; j<pointsInShower.size(); j++){ |
249 | double cell_E = pointsInShower[j]->E(); |
250 | double cell_r = pointsInShower[j]->r(); |
251 | double cell_phi = pointsInShower[j]->phi(); |
252 | E += cell_E; |
253 | x += cell_E*cell_r*cos(cell_phi); |
254 | sig_x += cell_E*cell_r*cos(cell_phi)*cell_r*cos(cell_phi); |
255 | y += cell_E*cell_r*sin(cell_phi); |
256 | sig_y += cell_E*cell_r*sin(cell_phi)*cell_r*sin(cell_phi); |
257 | |
258 | |
259 | double z_wt = 1/(pointsInShower[j]->sigZ()*pointsInShower[j]->sigZ()); |
260 | double cell_z = pointsInShower[j]->z(); |
261 | double cell_t = pointsInShower[j]->t(); |
262 | |
263 | sum_z_wt += z_wt; |
264 | z += z_wt*cell_z; |
265 | sig_z += z_wt*cell_z*cell_z; |
266 | |
267 | t += cell_E*cell_t; |
268 | sig_t += cell_E*cell_t*cell_t; |
269 | N_cell++; |
270 | |
271 | shower->AddAssociatedObject(pointsInShower[j]); |
272 | } |
273 | |
274 | x /= E; |
275 | sig_x /= E; |
276 | sig_x = sqrt(sig_x - x*x)/sqrt(N_cell); |
277 | y /= E; |
278 | sig_y /= E; |
279 | sig_y = sqrt(sig_y - y*y)/sqrt(N_cell); |
| Value stored to 'sig_y' is never read |
280 | z /= sum_z_wt; |
281 | sig_z /= sum_z_wt; |
282 | sig_z = sqrt(sig_z - z*z)/sqrt(N_cell); |
283 | t /= E; |
284 | sig_t /= E; |
285 | sig_t = sqrt(sig_t - t*t)/sqrt(N_cell); |
286 | |
287 | shower->id = id++; |
288 | shower->E_raw = E; |
289 | shower->x = x; |
290 | shower->y = y; |
291 | shower->z = z + m_z_target_center; |
292 | shower->t = t; |
293 | shower->N_cell = N_cell; |
294 | |
295 | |
296 | |
297 | |
298 | |
299 | |
300 | |
301 | |
302 | |
303 | |
304 | |
305 | |
306 | |
307 | float r = sqrt( shower->x * shower->x + shower->y * shower->y ); |
308 | |
309 | float zEntry = ( shower->z - m_z_target_center ) * ( bcalGeom->GetBCAL_inner_rad() / r ); |
310 | |
311 | float scale = m_scaleZ_p0 + m_scaleZ_p1*zEntry + |
312 | m_scaleZ_p2*(zEntry*zEntry) + m_scaleZ_p3*(zEntry*zEntry*zEntry); |
313 | float nonlin = m_nonlinZ_p0 + m_nonlinZ_p1*zEntry + |
314 | m_nonlinZ_p2*(zEntry*zEntry) + m_nonlinZ_p3*(zEntry*zEntry*zEntry); |
315 | |
316 | shower->E = pow( (shower->E_raw ) / scale, 1 / ( 1 + nonlin ) ); |
317 | |
318 | |
319 | |
320 | |
321 | |
322 | |
323 | |
324 | _data.push_back(shower); |
325 | } |
326 | |
327 | return NOERROR; |
328 | } |
329 | |
330 | |
331 | |
332 | |
333 | |
334 | void DBCALShower_factory_KLOE::FindPointsInShower(int indx, JEventLoop *loop, vector<const DBCALPoint*> &pointsInShower) |
335 | { |
336 | |
337 | |
338 | |
339 | |
340 | |
341 | |
342 | |
343 | |
344 | |
345 | |
346 | |
347 | |
348 | |
349 | |
350 | |
351 | |
352 | vector<const DBCALPoint*> points; |
353 | loop->Get(points); |
354 | |
355 | int start_indx = indx; |
356 | do{ |
357 | int module = narr[1][indx]; |
358 | int layer = narr[2][indx]; |
359 | int sector = narr[3][indx]; |
360 | |
361 | |
362 | for(unsigned int i=0; i<points.size(); i++){ |
363 | if(points[i]->module() !=module)continue; |
364 | if(points[i]->layer() !=layer)continue; |
365 | if(points[i]->sector() !=sector)continue; |
366 | pointsInShower.push_back(points[i]); |
367 | } |
368 | |
369 | |
370 | indx = next[indx]; |
371 | }while(indx != start_indx); |
372 | |
373 | } |
374 | |
375 | |
376 | |
377 | |
378 | |
379 | void DBCALShower_factory_KLOE::CellRecon(JEventLoop *loop) |
380 | { |
381 | |
382 | |
383 | |
384 | |
385 | |
386 | |
387 | |
388 | |
389 | |
390 | |
391 | |
392 | |
393 | |
394 | |
395 | |
396 | |
397 | |
398 | |
399 | |
400 | |
401 | memset( ecel_a, 0, modulemax_bcal48 * layermax_bcal10 * |
402 | colmax_bcal4 * sizeof( float ) ); |
403 | memset( tcel_a, 0, modulemax_bcal48 * layermax_bcal10 * |
404 | colmax_bcal4 * sizeof( float ) ); |
405 | memset( ecel_b, 0, modulemax_bcal48 * layermax_bcal10 * |
406 | colmax_bcal4 * sizeof( float ) ); |
407 | memset( tcel_b, 0, modulemax_bcal48 * layermax_bcal10 * |
408 | colmax_bcal4 * sizeof( float ) ); |
409 | |
410 | |
411 | |
412 | |
413 | |
414 | |
415 | |
416 | vector<const DBCALPoint*> points; |
417 | loop->Get(points); |
418 | if(points.size() <=0) return; |
419 | |
420 | for (vector<const DBCALPoint*>::const_iterator point_iter = points.begin(); |
421 | point_iter != points.end(); |
422 | ++point_iter) { |
423 | const DBCALPoint &point = **point_iter; |
424 | int module = point.module(); |
425 | int layer = point.layer(); |
426 | int sector = point.sector(); |
427 | double r = point.r(); |
428 | double phi = point.phi(); |
429 | double x = r*cos(phi); |
430 | double y = r*sin(phi); |
431 | |
432 | xcel[module-1][layer-1][sector-1] = x; |
433 | ycel[module-1][layer-1][sector-1] = y; |
434 | |
435 | |
436 | zcel[module-1][layer-1][sector-1] = point.z()+m_z_target_center-zOffset; |
437 | tcel[module-1][layer-1][sector-1] = point.t(); |
438 | ecel[module-1][layer-1][sector-1] = point.E(); |
439 | |
440 | |
441 | |
442 | double EUp=0,EDown=0,tUp=0,tDown=0; |
443 | vector<const DBCALUnifiedHit*> assoc_hits; |
444 | point.Get(assoc_hits); |
445 | for (unsigned int i=0; i<assoc_hits.size(); i++) { |
446 | if (assoc_hits[i]->end == DBCALGeometry::kUpstream) { |
447 | EUp = assoc_hits[i]->E; |
448 | tUp = assoc_hits[i]->t; |
449 | } |
450 | if (assoc_hits[i]->end == DBCALGeometry::kDownstream) { |
451 | EDown = assoc_hits[i]->E; |
452 | tDown = assoc_hits[i]->t; |
453 | } |
454 | |
455 | } |
456 | |
457 | ecel_a[module-1][layer-1][sector-1] = EUp; |
458 | |
459 | tcel_a[module-1][layer-1][sector-1] = tUp; |
460 | tcell_anor[module-1][layer-1][sector-1] = tUp; |
461 | |
462 | ecel_b[module-1][layer-1][sector-1] = EDown; |
463 | tcel_b[module-1][layer-1][sector-1] = tDown; |
464 | tcell_bnor[module-1][layer-1][sector-1] = tDown; |
465 | } |
466 | } |
467 | |
468 | |
469 | |
470 | |
471 | |
472 | void DBCALShower_factory_KLOE::CeleToArray(void) |
473 | { |
474 | |
475 | |
476 | |
477 | |
478 | |
479 | |
480 | |
481 | |
482 | |
483 | |
484 | |
485 | |
486 | |
487 | |
488 | |
489 | |
490 | |
491 | |
492 | |
493 | |
494 | |
495 | |
496 | |
497 | |
498 | |
499 | |
500 | |
501 | |
502 | |
503 | |
504 | |
505 | |
506 | |
507 | |
508 | |
509 | |
510 | celtot=0; |
511 | |
512 | for (int k = 0; k < modulemax_bcal48; k++){ |
513 | for (int i = 0; i < layermax_bcal10; i++){ |
514 | for (int j = 0; j < colmax_bcal4; j++){ |
515 | |
516 | float ea = ecel_a[k][i][j]; |
517 | float eb = ecel_b[k][i][j]; |
518 | float ta = tcel_a[k][i][j]; |
519 | float tb = tcel_b[k][i][j]; |
520 | |
521 | if( (min(ea,eb)>ethr_cell) & (fabs(ta-tb)<35.) & (ta!=0.) & (tb!=0.)) { |
522 | celtot=celtot+1; |
523 | } else { |
524 | continue; |
525 | } |
526 | |
527 | |
528 | if(celtot>cellmax_bcal48*10*4) { |
529 | break; |
530 | } |
531 | |
532 | narr[1][celtot]=k+1; |
533 | narr[2][celtot]=i+1; |
534 | narr[3][celtot]=j+1; |
535 | |
536 | |
537 | |
538 | celdata[1][celtot]=ea/0.145; |
539 | celdata[2][celtot]=eb/0.145; |
540 | |
541 | nclus[celtot] = celtot; |
542 | next[celtot] = celtot; |
543 | |
544 | e_cel[celtot] = ecel[k][i][j]; |
545 | x_cel[celtot] = xcel[k][i][j]; |
546 | y_cel[celtot] = ycel[k][i][j]; |
547 | z_cel[celtot] = zcel[k][i][j]; |
548 | t_cel[celtot] = tcel[k][i][j]; |
549 | |
550 | ta_cel[celtot]=tcell_anor[k][i][j]; |
551 | tb_cel[celtot]=tcell_bnor[k][i][j]; |
552 | } |
553 | } |
554 | } |
555 | } |
556 | |
557 | |
558 | |
559 | |
560 | |
561 | |
562 | void DBCALShower_factory_KLOE::PreCluster(JEventLoop *loop) |
563 | { |
564 | |
565 | |
566 | |
567 | |
568 | |
569 | |
570 | |
571 | |
572 | int k=1; |
573 | |
574 | |
575 | |
576 | |
577 | |
578 | |
579 | |
580 | |
581 | int modmin = 0; |
582 | int modmax = bcalGeom->GetBCAL_Nmodules(); |
583 | |
584 | |
585 | |
586 | int rowmax1= bcalGeom->GetBCAL_NInnerLayers(); |
587 | int rowmin2= rowmax1+1; |
588 | |
589 | int colmax1=bcalGeom->GetBCAL_NInnerSectors(); |
590 | int colmax2=bcalGeom->GetBCAL_NOuterSectors(); |
591 | |
592 | float r_middle= bcalGeom->GetBCAL_middle_rad(); |
593 | |
594 | |
595 | float thick_inner=bcalGeom->rSize(bcalGeom->cellId(1,bcalGeom->GetBCAL_NInnerLayers(),1)); |
596 | |
597 | float thick_outer=bcalGeom->rSize(bcalGeom->cellId(1,bcalGeom->GetBCAL_NInnerLayers()+1,1)); |
598 | |
599 | |
600 | float dis_in_out=bcalGeom->r(bcalGeom->cellId(1,bcalGeom->GetBCAL_NInnerLayers()+1,1))-bcalGeom->r(bcalGeom->cellId(1,bcalGeom->GetBCAL_NInnerLayers(),1)); |
601 | |
602 | float degree_permodule=360.0/(modmax-modmin); |
603 | float half_degree_permodule=degree_permodule/2.0; |
604 | |
605 | |
606 | float width_1=2.0*(r_middle-thick_inner/2.0)* |
607 | sin(half_degree_permodule*3.141593/180)/colmax1; |
608 | |
609 | float width_2=2.0*(r_middle+thick_outer/2.0)* |
610 | sin(half_degree_permodule*3.141593/180)/colmax2; |
611 | |
612 | |
613 | |
614 | float disthres=width_2*1.5-width_1*0.5+0.0001; |
615 | |
616 | for (int i = 1; i < (celtot+1); i++){ |
617 | |
618 | int maxnn=0; |
619 | float emin=0.; |
620 | |
621 | |
622 | for (int j = 1; j < (celtot+1); j++){ |
623 | if ( (j!=i) & (nclus[j]!=nclus[i]) & (e_cel[j]>emin)) { |
624 | |
625 | int k1= narr[1][i]; |
626 | int k2= narr[1][j]; |
627 | int i1= narr[2][i]; |
628 | int i2= narr[2][j]; |
629 | |
630 | |
631 | int modiff = k1-k2; |
632 | int amodif = abs(modiff); |
633 | |
634 | |
635 | if ( (abs(i1-i2)<=k) & ((amodif<=1) || (amodif==47)) ) { |
636 | |
637 | int j1= narr[3][i]; |
638 | int j2= narr[3][j]; |
639 | |
640 | if(amodif==0) { |
641 | |
642 | if ( (i1<=rowmax1) & (i2<=rowmax1) & (abs(j2-j1)<=k) ) { |
643 | emin=e_cel[j]; |
644 | maxnn=j; |
645 | } |
646 | |
647 | |
648 | |
649 | if ( (i1>=rowmin2) & (i2>=rowmin2) & (abs(j2-j1)<=k) ) { |
650 | emin=e_cel[j]; |
651 | maxnn=j; |
652 | } |
653 | } |
654 | |
655 | if(amodif>0) { |
656 | if( (modiff==1) || (modiff==-47) ) { |
657 | if ( (i1<=rowmax1) & (i2<=rowmax1) ){ |
658 | if(abs((j1+colmax1)-j2)<=k){ |
659 | emin=e_cel[j]; |
660 | maxnn=j; |
661 | } |
662 | } |
663 | |
664 | if ( (i1>=rowmin2) & (i2>=rowmin2) ) { |
665 | if(abs((j1+colmax2)-j2)<=k){ |
666 | emin=e_cel[j]; |
667 | maxnn=j; |
668 | } |
669 | } |
670 | } |
671 | |
672 | if ( (modiff==-1) || (modiff==47) ) { |
673 | |
674 | if ( (i1<=rowmax1) & (i2<=rowmax1) ){ |
675 | if(abs((j2+colmax1)-j1)<=k){ |
676 | emin=e_cel[j]; |
677 | maxnn=j; |
678 | } |
679 | } |
680 | |
681 | if ( (i1>=rowmin2) & (i2>=rowmin2) ){ |
682 | if(abs((j2+colmax2)-j1)<=k){ |
683 | emin=e_cel[j]; |
684 | maxnn=j; |
685 | } |
686 | } |
687 | } |
688 | } |
689 | |
690 | |
691 | |
692 | |
693 | if( ( (i1 == rowmax1) & (i2 == rowmin2) ) || |
694 | ( (i1 == rowmin2) & (i2 == rowmax1) ) ) { |
695 | |
696 | float delta_xx=xx[k1-1][i1-1][j1-1]-xx[k2-1][i2-1][j2-1]; |
697 | float delta_yy=yy[k1-1][i1-1][j1-1]-yy[k2-1][i2-1][j2-1]; |
698 | |
699 | |
700 | float dis = sqrt( delta_xx * delta_xx + delta_yy * delta_yy ); |
701 | |
702 | dis = sqrt( dis*dis - dis_in_out * dis_in_out ); |
703 | |
704 | if( dis < disthres ){ |
705 | emin = e_cel[j]; |
706 | maxnn = j; |
707 | } |
708 | } |
709 | } |
710 | } |
711 | } |
712 | |
713 | if(maxnn>0){ |
714 | |
715 | Connect(maxnn,i); |
716 | } |
717 | } |
718 | } |
719 | |
720 | |
721 | |
722 | |
723 | |
724 | |
725 | void DBCALShower_factory_KLOE::Connect(int n,int m) |
726 | { |
727 | |
728 | |
729 | |
730 | |
731 | |
732 | |
733 | |
734 | |
735 | |
736 | |
737 | |
738 | |
739 | |
740 | |
741 | |
742 | |
743 | |
744 | |
745 | |
746 | |
747 | |
748 | |
749 | |
750 | |
751 | |
752 | |
753 | |
754 | |
755 | |
756 | |
757 | |
758 | |
759 | |
760 | if(nclus[n]!=nclus[m]){ |
761 | int j=m; |
762 | nclus[j]=nclus[n]; |
763 | while(next[j]!=m){ |
764 | j=next[j]; |
765 | nclus[j]=nclus[n]; |
766 | } |
767 | next[j]=next[n]; |
768 | next[n]=m; |
769 | } |
770 | } |
771 | |
772 | |
773 | |
774 | |
775 | |
776 | void DBCALShower_factory_KLOE::ClusNorm(void) |
777 | { |
778 | |
779 | memset( e_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
780 | memset( x_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
781 | memset( y_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
782 | memset( z_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
783 | memset( t_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
784 | memset( ea_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
785 | memset( eb_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
786 | memset( ta_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
787 | memset( tb_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
788 | memset( tsqr_a, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
789 | memset( tsqr_b, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
790 | memset( trms_a, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
791 | memset( trms_b, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
792 | memset( e2_a, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
793 | memset( e2_b, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
794 | memset( clspoi, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
795 | memset( ncltot, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
796 | memset( ntopol, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
797 | |
798 | |
799 | |
800 | |
801 | |
802 | |
803 | |
804 | |
805 | |
806 | |
807 | clstot=0; |
808 | |
809 | for (int ix = 1; ix < (celtot+1); ix++){ |
810 | |
811 | |
812 | |
813 | |
814 | |
815 | int n=nclus[ix]; |
816 | int j=0; |
817 | |
818 | for (int i = 1; i < (clstot+1); i++){ |
819 | if(n==clspoi[i]) j=i; |
820 | } |
821 | |
822 | if(j==0) { |
823 | clstot=clstot+1; |
824 | clspoi[clstot]=n; |
825 | } |
826 | |
827 | |
828 | |
829 | |
830 | if(e_cel[ix]<0.000000001)continue; |
831 | |
832 | x_cls[n]=(e_cls[n]*x_cls[n]+e_cel[ix]*x_cel[ix]) |
833 | /(e_cls[n]+e_cel[ix]); |
834 | |
835 | y_cls[n]=(e_cls[n]*y_cls[n]+e_cel[ix]*y_cel[ix]) |
836 | /(e_cls[n]+e_cel[ix]); |
837 | |
838 | z_cls[n]=(e_cls[n]*z_cls[n]+e_cel[ix]*z_cel[ix]) |
839 | /(e_cls[n]+e_cel[ix]); |
840 | |
841 | t_cls[n]=(e_cls[n]*t_cls[n]+e_cel[ix]*t_cel[ix]) |
842 | /(e_cls[n]+e_cel[ix]); |
843 | |
844 | e_cls[n]=e_cls[n]+e_cel[ix]; |
845 | |
846 | |
847 | |
848 | |
849 | |
850 | |
851 | ta_cls[n]=(ea_cls[n]*ta_cls[n]+celdata[1][ix]*ta_cel[ix]) |
852 | /(ea_cls[n]+celdata[1][ix]); |
853 | tsqr_a[n]=(ea_cls[n]*tsqr_a[n]+celdata[1][ix]*ta_cel[ix]* |
854 | ta_cel[ix])/(ea_cls[n]+celdata[1][ix]); |
855 | ea_cls[n]=ea_cls[n]+celdata[1][ix]; |
856 | e2_a[n]=e2_a[n]+celdata[1][ix]*celdata[1][ix]; |
857 | tb_cls[n]=(eb_cls[n]*tb_cls[n]+celdata[2][ix]*tb_cel[ix])/ |
858 | (eb_cls[n]+celdata[2][ix]); |
859 | tsqr_b[n]=(eb_cls[n]*tsqr_b[n]+celdata[2][ix]*tb_cel[ix]* |
860 | tb_cel[ix])/(eb_cls[n]+celdata[2][ix]); |
861 | eb_cls[n]=eb_cls[n]+celdata[2][ix]; |
862 | e2_b[n]=e2_b[n]+celdata[2][ix]*celdata[2][ix]; |
863 | |
864 | |
865 | |
866 | |
867 | |
868 | ncltot[n]++; |
869 | |
870 | if( narr[1][n] != narr[1][ix] || narr[2][n] != narr[2][ix] ) |
871 | ntopol[n]++; |
872 | |
873 | |
874 | |
875 | } |
876 | |
877 | for (int n = 1; n < (clstot+1); n++){ |
878 | |
879 | int ix=clspoi[n]; |
880 | if( ncltot[ix] > 1) { |
881 | |
882 | float effnum = ea_cls[ix] * ea_cls[ix] / e2_a[ix]; |
883 | trms_a[ix] = effnum / ( effnum - 1 ) * |
884 | ( tsqr_a[ix] - ta_cls[ix] * ta_cls[ix] ); |
885 | |
886 | effnum = eb_cls[ix] * eb_cls[ix] / e2_b[ix]; |
887 | trms_b[ix] = effnum / ( effnum - 1 ) * |
888 | ( tsqr_b[ix] - tb_cls[ix] * tb_cls[ix] ); |
889 | |
890 | if( trms_a[ix] <= 0.0 ) trms_a[ix] = 0.; |
891 | if( trms_b[ix] <= 0.0 ) trms_b[ix] = 0.; |
892 | trms_a[ix] = sqrt( trms_a[ix] ); |
893 | trms_b[ix] = sqrt( trms_b[ix] ); |
894 | } |
895 | else { |
896 | trms_a[ix] = 0.; |
897 | trms_b[ix] = 0.; |
898 | } |
899 | } |
900 | } |
901 | |
902 | |
903 | |
904 | |
905 | void DBCALShower_factory_KLOE::ClusAnalysis() |
906 | { |
907 | |
908 | |
909 | bool newClust = false; |
910 | |
911 | |
912 | |
913 | |
914 | |
915 | for (int i = 0; i < 2; i++){ |
916 | for (int j = 1; j < (clstot+1); j++){ |
917 | int ix=clspoi[j]; |
918 | if(e_cls[ix]>0.0){ |
919 | float dist=sqrt(trms_a[ix]*trms_a[ix]+trms_b[ix]*trms_b[ix]); |
920 | if(dist>BREAK_THRESH_TRMS) { |
921 | Clus_Break(ix); |
922 | newClust = true; |
923 | } |
924 | } |
925 | } |
926 | |
927 | if( newClust ){ |
928 | |
929 | ClusNorm(); |
930 | newClust = false; |
931 | } |
932 | } |
933 | |
934 | |
935 | |
936 | |
937 | |
938 | int icls[3]; |
939 | for (int i = 1; i < clstot; i++){ |
940 | icls[1]=0; |
941 | icls[2]=0; |
942 | for (int j = (i+1); j < (clstot+1); j++){ |
943 | |
944 | int ix=clspoi[i]; |
945 | int iy=clspoi[j]; |
946 | |
947 | |
948 | if ( (e_cls[ix]>0.0) & (e_cls[iy]>0.0) ) { |
949 | |
950 | float delta_x=x_cls[ix]-x_cls[iy]; |
951 | float delta_y=y_cls[ix]-y_cls[iy]; |
952 | float delta_z=z_cls[ix]-z_cls[iy]; |
953 | float dist=sqrt(delta_x*delta_x+delta_y*delta_y+delta_z*delta_z); |
954 | |
955 | float tdif=fabs(t_cls[ix]-t_cls[iy]); |
956 | |
957 | |
958 | |
959 | |
960 | if ( (dist<MERGE_THRESH_DIST) & (tdif<MERGE_THRESH_TIME) ){ |
961 | float zdif=fabs(z_cls[ix]-z_cls[iy]); |
962 | float distran=sqrt(delta_x*delta_x+delta_y*delta_y); |
963 | |
964 | if ( (zdif<MERGE_THRESH_ZDIST) & (distran<MERGE_THRESH_XYDIST) ){ |
965 | if(e_cls[ix]>=e_cls[iy]) { |
966 | icls[1]=ix; |
967 | icls[2]=iy; |
968 | } |
969 | else { |
970 | icls[1]=iy; |
971 | icls[2]=ix; |
972 | } |
973 | } |
974 | } |
975 | |
976 | } |
977 | |
978 | if(min(icls[1],icls[2])>0){ |
979 | |
980 | Connect(icls[1],icls[2]); |
981 | newClust = true; |
982 | } |
983 | } |
984 | } |
985 | |
986 | if( newClust ){ |
987 | |
988 | ClusNorm(); |
989 | } |
990 | } |
991 | |
992 | |
993 | |
994 | |
995 | void DBCALShower_factory_KLOE::Clus_Break(int nclust) |
996 | { |
997 | int nseed[5],selnum,selcel[cellmax_bcal48*10*4+1]; |
998 | float tdif,tdif_a,tdif_b,tseed[5]; |
999 | |
1000 | |
1001 | for (int i =0; i < 5; i++){ |
1002 | nseed[i]=0; |
1003 | tseed[i]=0; |
1004 | } |
1005 | |
1006 | |
1007 | |
1008 | int n=nclust; |
1009 | tdif_a=ta_cel[n]-ta_cls[nclust]; |
1010 | tdif_b=tb_cel[n]-tb_cls[nclust]; |
1011 | |
1012 | |
1013 | |
1014 | if(tdif_a>0.0) { |
1015 | if(tdif_b>0){ |
1016 | selnum=1; |
1017 | } |
1018 | else { |
1019 | selnum=2; |
1020 | } |
1021 | } |
1022 | |
1023 | else { |
1024 | if(tdif_b>0.0) { |
1025 | selnum=3; |
1026 | } |
1027 | else { |
1028 | selnum=4; |
1029 | } |
1030 | } |
1031 | |
1032 | |
1033 | |
1034 | if(selnum>0) { |
1035 | float tdif=sqrt(tdif_a*tdif_a+tdif_b*tdif_b); |
1036 | if(tdif>tseed[selnum]){ |
1037 | nseed[selnum]=n; |
1038 | tseed[selnum]=tdif; |
1039 | } |
1040 | selcel[n]=selnum; |
1041 | } |
1042 | |
1043 | |
1044 | |
1045 | while(next[n]!=nclust) { |
1046 | n=next[n]; |
1047 | tdif_a=ta_cel[n]-ta_cls[nclust]; |
1048 | tdif_b=tb_cel[n]-tb_cls[nclust]; |
1049 | |
1050 | |
1051 | |
1052 | |
1053 | if(tdif_a>0.0) { |
1054 | |
1055 | if(tdif_b>0.0) { |
1056 | selnum=1; |
1057 | } |
1058 | else { |
1059 | selnum=2; |
1060 | } |
1061 | } |
1062 | |
1063 | else { |
1064 | if(tdif_b>0.0) { |
1065 | selnum=3; |
1066 | } |
1067 | else { |
1068 | selnum=4; |
1069 | } |
1070 | } |
1071 | |
1072 | |
1073 | |
1074 | |
1075 | |
1076 | |
1077 | if(selnum>0){ |
1078 | tdif=sqrt(tdif_a*tdif_a+tdif_b*tdif_b); |
1079 | |
1080 | if(tdif>tseed[selnum]){ |
1081 | nseed[selnum]=n; |
1082 | tseed[selnum]=tdif; |
1083 | } |
1084 | |
1085 | selcel[n]=selnum; |
1086 | } |
1087 | |
1088 | |
1089 | |
1090 | } |
1091 | |
1092 | |
1093 | |
1094 | |
1095 | |
1096 | |
1097 | for (int i =1; i < 5; i++){ |
1098 | |
1099 | if(nseed[i]>0) { |
1100 | |
1101 | nclus[nseed[i]]=nseed[i]; |
1102 | next[nseed[i]]=nseed[i]; |
1103 | |
1104 | for (int j =1; j < (celtot+1); j++){ |
1105 | if ( (nclus[j]==nclust) & (j!=nseed[i]) ){ |
1106 | if(selcel[j]==i) { |
1107 | nclus[j]=j; |
1108 | next[j]=j; |
1109 | Connect(nseed[i],j); |
1110 | } |
1111 | } |
1112 | } |
1113 | } |
1114 | } |
1115 | } |
1116 | |
1117 | |
1118 | |
1119 | |
1120 | |
1121 | |
1122 | |
1123 | |
1124 | |
1125 | void DBCALShower_factory_KLOE::Trakfit( void ) |
1126 | { |
1127 | |
1128 | float emin=0.0001; |
1129 | |
1130 | memset( clslyr, 0, ( clsmax_bcal48*10*4 + 1 ) * |
1131 | ( layermax_bcal10 + 1 ) * 6 * sizeof( float ) ); |
1132 | |
1133 | memset( apx, 0, ( clsmax_bcal48*10*4 + 1 ) * 4 * sizeof( float ) ); |
1134 | memset( eapx, 0, ( clsmax_bcal48*10*4 + 1 ) * 4 * sizeof( float ) ); |
1135 | memset( ctrk, 0, ( clsmax_bcal48*10*4 + 1 ) * 4 * sizeof( float ) ); |
1136 | memset( ectrk, 0, ( clsmax_bcal48*10*4 + 1 ) * 4 * sizeof( float ) ); |
1137 | |
1138 | for (int ix = 1; ix < (celtot+1); ix++){ |
1139 | |
1140 | int n = nclus[ix]; |
1141 | |
1142 | |
1143 | |
1144 | |
1145 | |
1146 | int lyr = narr[2][ix]; |
1147 | |
1148 | |
1149 | clslyr[xlyr][lyr][n]=(clslyr[elyr][lyr][n]*clslyr[xlyr][lyr][n] |
1150 | +e_cel[ix]*x_cel[ix])/(clslyr[elyr][lyr][n]+e_cel[ix]); |
1151 | |
1152 | clslyr[ylyr][lyr][n]=(clslyr[elyr][lyr][n]*clslyr[ylyr][lyr][n] |
1153 | +e_cel[ix]*y_cel[ix])/(clslyr[elyr][lyr][n]+e_cel[ix]); |
1154 | |
1155 | |
1156 | clslyr[zlyr][lyr][n]=(clslyr[elyr][lyr][n]*clslyr[zlyr][lyr][n] |
1157 | +e_cel[ix]*z_cel[ix])/(clslyr[elyr][lyr][n]+e_cel[ix]); |
1158 | |
1159 | |
1160 | clslyr[tlyr][lyr][n]=(clslyr[elyr][lyr][n]*clslyr[tlyr][lyr][n] |
1161 | +e_cel[ix]*t_cel[ix])/(clslyr[elyr][lyr][n]+e_cel[ix]); |
1162 | |
1163 | |
1164 | clslyr[elyr][lyr][n]=clslyr[elyr][lyr][n]+e_cel[ix]; |
1165 | |
1166 | |
1167 | } |
1168 | |
1169 | memset( nlrtot, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) ); |
1170 | |
1171 | for (int n = 1; n < ( clstot + 1 ); n++){ |
1172 | |
1173 | int ix=clspoi[n]; |
1174 | |
1175 | for (int i = 1; i < (layermax_bcal10+1); i++){ |
1176 | |
1177 | if( clslyr[elyr][i][ix] > 0.0 ) nlrtot[ix]++; |
1178 | } |
1179 | |
1180 | for (int i = 0; i < ( layermax_bcal10 + 1 ); i++){ |
1181 | |
1182 | x[i]=0.0; |
1183 | y[i]=0.0; |
1184 | z[i]=0.0; |
1185 | e[i]=0.0; |
1186 | sigx[i]=0.0; |
1187 | sigy[i]=0.0; |
1188 | sigz[i]=0.0; |
1189 | } |
1190 | |
1191 | int nltot=0; |
1192 | |
1193 | for (int il = 1; il < (layermax_bcal10+1); il++){ |
1194 | |
1195 | if(clslyr[elyr][il][ix]>emin) { |
1196 | |
1197 | nltot=nltot+1; |
1198 | x[nltot]= clslyr[xlyr][il][ix]; |
1199 | y[nltot]= clslyr[ylyr][il][ix]; |
1200 | z[nltot]= clslyr[zlyr][il][ix]; |
1201 | e[nltot]= clslyr[elyr][il][ix]; |
1202 | |
1203 | sigy[nltot] = 1.0/e[nltot]; |
1204 | sigx[nltot] = 1.0/e[nltot]; |
1205 | sigz[nltot] = 1.0/sqrt(e[nltot]); |
1206 | } |
1207 | } |
1208 | |
1209 | |
1210 | |
1211 | |
1212 | |
1213 | |
1214 | |
1215 | |
1216 | sigx[1]=0.5; |
1217 | sigx[2]=0.5; |
1218 | sigx[3]=0.5; |
1219 | sigx[4]=0.5; |
1220 | sigx[5]=0.5; |
1221 | sigx[6]=0.8; |
1222 | sigx[7]=0.9; |
1223 | sigx[8]=1.2; |
1224 | sigx[9]=1.3; |
1225 | |
1226 | sigy[1]=0.5; |
1227 | sigy[2]=0.5; |
1228 | sigy[3]=0.5; |
1229 | sigy[4]=0.5; |
1230 | sigy[5]=0.5; |
1231 | sigy[6]=0.8; |
1232 | sigy[7]=0.9; |
1233 | sigy[8]=1.2; |
1234 | sigy[9]=1.3; |
1235 | |
1236 | |
1237 | sigz[1]=0.5; |
1238 | sigz[2]=0.5; |
1239 | sigz[3]=0.5; |
1240 | sigz[4]=0.5; |
1241 | sigz[5]=0.5; |
1242 | sigz[6]=0.8; |
1243 | sigz[7]=0.9; |
1244 | sigz[8]=1.2; |
1245 | sigz[9]=1.3; |
1246 | |
1247 | if( nltot > 1 ){ |
1248 | |
1249 | Fit_ls(); |
1250 | |
1251 | for (int i = 1; i < 4; i++){ |
1252 | |
1253 | ctrk[i][ix]=ctrk_ix[i]; |
1254 | ectrk[i][ix]=ectrk_ix[i]; |
1255 | apx[i][ix]=apx_ix[i]; |
1256 | eapx[i][ix]=eapx_ix[i]; |
1257 | } |
1258 | } |
1259 | else{ |
1260 | |
1261 | |
1262 | |
1263 | apx[1][ix] = x[1]; |
1264 | apx[2][ix] = y[1]; |
1265 | apx[3][ix] = z[1]; |
1266 | eapx[1][ix] = sigx[1]; |
1267 | eapx[2][ix] = sigy[1]; |
1268 | eapx[3][ix] = sigz[1]; |
1269 | ectrk[1][ix] = 0.0; |
1270 | ectrk[2][ix] = 0.0; |
1271 | ectrk[3][ix] = 0.0; |
1272 | ctrk[1][ix] = 999.0; |
1273 | ctrk[2][ix] = 999.0; |
1274 | ctrk[3][ix] = 999.0; |
1275 | } |
1276 | } |
1277 | } |
1278 | |
1279 | |
1280 | |
1281 | |
1282 | void DBCALShower_factory_KLOE::Fit_ls() |
1283 | { |
1284 | float a,b,c; |
1285 | float d,e,f,chi2,q,norm; |
1286 | float siga,sigb,sigc,sigd,sige,sigf; |
1287 | float sigb2,sigd2,sigf2; |
1288 | |
1289 | |
1290 | Linefit(1,1,a,b,siga,sigb,chi2,q); |
1291 | |
1292 | Linefit(2,1,c,d,sigc,sigd,chi2,q); |
1293 | |
1294 | Linefit(3,1,e,f,sige,sigf,chi2,q); |
1295 | sigb2=sigb*sigb; |
1296 | sigd2=sigd*sigd; |
1297 | sigf2=sigf*sigf; |
1298 | |
1299 | apx_ix[1]=a; |
1300 | apx_ix[2]=c; |
1301 | apx_ix[3]=e; |
1302 | eapx_ix[1]=siga; |
1303 | eapx_ix[2]=sigc; |
1304 | eapx_ix[3]=sige; |
1305 | |
1306 | |
1307 | |
1308 | |
1309 | |
1310 | |
1311 | norm=sqrt(b*b+d*d+f*f); |
1312 | |
1313 | ctrk_ix[1]=b/norm; |
1314 | ctrk_ix[2]=d/norm; |
1315 | ctrk_ix[3]=f/norm; |
1316 | |
1317 | float norm3=norm*norm*norm; |
1318 | |
1319 | ectrk_ix[1]=sqrt((d*d+f*f)*(d*d+f*f)*sigb2+b*b*d*d*sigd2+b*b*f*f*sigf2)/norm3; |
1320 | ectrk_ix[2]=sqrt((b*b+f*f)*(b*b+f*f)*sigd2+d*d*b*b*sigb2+d*d*f*f*sigf2)/norm3; |
1321 | ectrk_ix[3]=sqrt((b*b+d*d)*(b*b+d*d)*sigf2+f*f*b*b*sigb2+f*f*d*d*sigd2)/norm3; |
1322 | |
1323 | return; |
1324 | } |
1325 | |
1326 | |
1327 | |
1328 | |
1329 | |
1330 | void DBCALShower_factory_KLOE::Linefit(int ixyz,int mwt,float &a, |
1331 | float &b,float &siga,float &sigb,float &chi2,float &q) |
1332 | { |
1333 | |
1334 | |
1335 | |
1336 | |
1337 | |
1338 | |
1339 | float sig[layermax_bcal10+1],etemp; |
1340 | float xtemp[layermax_bcal10+1],ytemp[layermax_bcal10+1]; |
1341 | |
1342 | |
1343 | |
1344 | |
1345 | |
1346 | |
1347 | |
1348 | |
1349 | |
1350 | |
1351 | |
1352 | |
1353 | float sigdat,ss,st2,sx,sxoss,sy,t,wt; |
1354 | sx=0.0; |
1355 | sy=0.0; |
1356 | st2=0.0; |
1357 | b=0.0; |
1358 | |
1359 | int ndata=0; |
1360 | |
1361 | |
1362 | |
1363 | if(ixyz==1) { |
1364 | for (int i = 1; i < (layermax_bcal10+1); i++){ |
1365 | xtemp[i]=rt[i]; |
1366 | ytemp[i]=x[i]; |
1367 | sig[i]=sigx[i]; |
1368 | etemp=e[i]; |
1369 | if(etemp>0.0001)ndata=ndata+1; |
1370 | |
1371 | } |
1372 | } |
1373 | else if(ixyz==2) { |
1374 | for (int i = 1; i < (layermax_bcal10+1); i++){ |
1375 | xtemp[i]=rt[i]; |
1376 | ytemp[i]=y[i]; |
1377 | sig[i]=sigy[i]; |
1378 | etemp=e[i]; |
1379 | if(etemp>0.000001)ndata=ndata+1; |
1380 | } |
1381 | } |
1382 | else if(ixyz==3) { |
1383 | for (unsigned int i = 1; i < (layermax_bcal10+1); i++){ |
1384 | xtemp[i]=rt[i]; |
1385 | ytemp[i]=z[i]; |
1386 | sig[i]=sigz[i]; |
1387 | etemp=e[i]; |
1388 | if(etemp>0.000001)ndata=ndata+1; |
1389 | } |
1390 | } |
1391 | |
1392 | if(mwt!=0) { |
1393 | ss=0.0; |
1394 | for (int i = 1; i < (ndata+1); i++){ |
1395 | wt=1.0/(sig[i]*sig[i]); |
1396 | ss=ss+wt; |
1397 | sx=sx+xtemp[i]*wt; |
1398 | sy=sy+ytemp[i]*wt; |
1399 | } |
1400 | } |
1401 | |
1402 | else { |
1403 | for (int i = 1; i < (ndata+1); i++){ |
1404 | sx=sx+xtemp[i]; |
1405 | sy=sy+ytemp[i]; |
1406 | } |
1407 | ss=float(ndata); |
1408 | } |
1409 | |
1410 | sxoss=sx/ss; |
1411 | |
1412 | if(mwt!=0) { |
1413 | for (int i = 1; i < (ndata+1); i++){ |
1414 | t=(xtemp[i]-sxoss)/sig[i]; |
1415 | st2=st2+t*t; |
1416 | b=b+t*ytemp[i]/sig[i]; |
1417 | } |
1418 | |
1419 | } |
1420 | |
1421 | else { |
1422 | |
1423 | for (int i = 1; i < (ndata+1); i++){ |
1424 | t=xtemp[i]-sxoss; |
1425 | st2=st2+t*t; |
1426 | b=b+t*ytemp[i]; |
1427 | } |
1428 | } |
1429 | |
1430 | b=b/st2; |
1431 | a=(sy-sx*b)/ss; |
1432 | siga=sqrt((1.0+sx*sx/(ss*st2))/ss); |
1433 | sigb=sqrt(1.0/st2); |
1434 | chi2=0.0; |
1435 | |
1436 | if(mwt==0) { |
1437 | for (int i = 1; i < (ndata+1); i++){ |
1438 | chi2=chi2+(ytemp[i]-a-b*xtemp[i])*(ytemp[i]-a-b*xtemp[i]); |
1439 | } |
1440 | q=1.0; |
1441 | sigdat=sqrt(chi2/(ndata-2)); |
1442 | siga=siga*sigdat; |
1443 | sigb=sigb*sigdat; |
1444 | } |
1445 | else { |
1446 | for (int i = 1; i < (ndata+1); i++){ |
1447 | chi2=chi2+((ytemp[i]-a-b*xtemp[i])/ |
1448 | sig[i])*((ytemp[i]-a-b*xtemp[i])/sig[i]); |
1449 | } |
1450 | q=Gammq(0.5*(ndata-2),0.5*chi2); |
1451 | } |
1452 | |
1453 | } |
1454 | |
1455 | |
1456 | |
1457 | |
1458 | |
1459 | |
1460 | float DBCALShower_factory_KLOE::Gammq(float a_gammq,float x_gammq) |
1461 | { |
1462 | |
1463 | |
1464 | |
1465 | |
1466 | float gammq; |
1467 | |
1468 | |
1469 | |
1470 | |
1471 | |
1472 | float gammcf=0,gamser; |
1473 | |
1474 | if(a_gammq==0.0) { |
1475 | gammq=999.0; |
1476 | return gammq; |
1477 | } |
1478 | |
1479 | |
1480 | |
1481 | if(x_gammq<0. || a_gammq<= 0.0) { |
1482 | |
1483 | return 999.0; |
1484 | } |
1485 | |
1486 | if(x_gammq<(a_gammq+1.)) { |
1487 | Gser(gamser,a_gammq,x_gammq); |
1488 | gammq=1.0-gamser; |
1489 | } |
1490 | else { |
1491 | Gcf(gammcf,a_gammq,x_gammq); |
1492 | gammq=gammcf; |
1493 | } |
1494 | return gammq; |
1495 | } |
1496 | |
1497 | |
1498 | |
1499 | |
1500 | |
1501 | void DBCALShower_factory_KLOE::Gser(float &gamser,float a_gser,float x_gser) |
1502 | { |
1503 | |
1504 | int itmax=100; |
1505 | float eps=3.0e-7; |
1506 | float gln; |
1507 | |
1508 | |
1509 | |
1510 | |
1511 | |
1512 | float ap, del,sum; |
1513 | |
1514 | gln=Gammln(a_gser); |
1515 | |
1516 | if(x_gser<=0.0) { |
1517 | if(x_gser<0.0) cout<<"x_gser<0 in gser"<<"\n"; |
1518 | gamser=0.0; |
1519 | return; |
1520 | } |
1521 | |
1522 | ap=a_gser; |
1523 | sum=1.0/a_gser; |
1524 | del=sum; |
1525 | |
1526 | |
1527 | for (int n = 1; n < (itmax+1); n++){ |
1528 | ap=ap+1.0; |
1529 | del=del*x_gser/ap; |
1530 | sum=sum+del; |
1531 | |
1532 | if(fabs(del)<fabs(sum)*eps) { |
1533 | gamser=sum*exp(-x_gser+a_gser*log(x_gser)-gln); |
1534 | return; |
1535 | } |
1536 | |
1537 | } |
1538 | |
1539 | |
1540 | return; |
1541 | } |
1542 | |
1543 | |
1544 | |
1545 | |
1546 | |
1547 | void DBCALShower_factory_KLOE::Gcf(float &gammcf,float a_gcf,float x_gcf) |
1548 | { |
1549 | |
1550 | |
1551 | int itmax=100; |
1552 | float eps=3.0e-7; |
1553 | float fpmin=1.0e-30; |
1554 | |
1555 | float gln; |
1556 | |
1557 | |
1558 | |
1559 | |
1560 | |
1561 | |
1562 | |
1563 | |
1564 | |
1565 | |
1566 | float an,b,c,d,del,h; |
1567 | |
1568 | gln=Gammln(a_gcf); |
1569 | b=x_gcf+1.0-a_gcf; |
1570 | c=1.0/fpmin; |
1571 | d=1.0/b; |
1572 | h=d; |
1573 | |
1574 | |
1575 | for (int i = 1; i < (itmax+1); i++){ |
1576 | an=-i*(i-a_gcf); |
1577 | b=b+2.0; |
1578 | d=an*d+b; |
1579 | if(fabs(d)<fpmin)d=fpmin; |
1580 | c=b+an/c; |
1581 | if(fabs(c)<fpmin)c=fpmin; |
1582 | d=1.0/d; |
1583 | del=d*c; |
1584 | h=h*del; |
1585 | if(fabs(del-1.0)<eps) { |
1586 | gammcf=exp(-x_gcf+a_gcf*log(x_gcf)-gln)*h; |
1587 | return; |
1588 | } |
1589 | |
1590 | return; |
1591 | } |
1592 | } |
1593 | |
1594 | |
1595 | |
1596 | |
1597 | float DBCALShower_factory_KLOE::Gammln(float xx_gln) |
1598 | { |
1599 | |
1600 | |
1601 | float ser,stp,tmp,x_gln,y_gln; |
1602 | float cof[7]; |
1603 | float gammln; |
1604 | |
1605 | |
1606 | |
1607 | |
1608 | |
1609 | |
1610 | |
1611 | |
1612 | stp=2.5066282746310005; |
1613 | cof[1]=76.18009172947146; |
1614 | cof[2]=-86.50532032941677; |
1615 | cof[3]=24.01409824083091; |
1616 | cof[4]=-1.231739572450155; |
1617 | cof[5]=.1208650973866179e-2; |
1618 | cof[6]=-.5395239384953e-5; |
1619 | |
1620 | x_gln=xx_gln; |
1621 | y_gln=x_gln; |
1622 | tmp=x_gln+5.5; |
1623 | tmp=(x_gln+0.5)*log(tmp)-tmp; |
1624 | |
1625 | ser=1.000000000190015; |
1626 | |
1627 | for (int j = 1; j < 7; j++){ |
1628 | y_gln=y_gln+1.0; |
1629 | ser=ser+cof[j]/y_gln; |
1630 | } |
1631 | |
1632 | |
1633 | gammln=tmp+log(stp*ser/x_gln); |
1634 | return gammln; |
1635 | } |
1636 | |