Bug Summary

File:libraries/BCAL/DBCALShower_factory_KLOE.cc
Location:line 1498, column 9
Description:Assigned value is garbage or undefined

Annotated Source Code

1// $Id: DBCALShower_factory_KLOE.cc 9234 2012-06-20 17:28:28Z wilevine $
2//
3// File: DBCALShower_factory_KLOE.cc
4// Created: Tue Jul 3 18:25:12 EDT 2007
5// Creator: Matthew Shepherd
6//
7
8#include <cassert>
9#include <math.h>
10#include <map>
11
12#include "BCAL/DBCALHit.h"
13#include "BCAL/DBCALGeometry.h"
14#include "BCAL/DBCALShower_factory_KLOE.h"
15
16#include "units.h"
17
18using namespace std;
19
20//------------------
21// DBCALShower_factory_KLOE
22//------------------
23DBCALShower_factory_KLOE::DBCALShower_factory_KLOE()
24{
25 // this should be lower than cut in mcsmear
26 ethr_cell=0.0001; // MIN ENERGY THRESD OF cell in GeV
27
28 CLUST_THRESH = 0.03; // MIN ENERGY THRESD OF CLUSTER IN GEV (make this match the value used in the other algorithm)
29
30 elyr = 1;
31 xlyr = 2;
32 ylyr = 3;
33 zlyr = 4;
34 tlyr = 5;
35
36 // these four parameters are used in merging clusters
37 MERGE_THRESH_DIST = 40.0; // CENTROID DISTANCE THRESHOLD
38 MERGE_THRESH_TIME = 2.5; // CENTROID TIME THRESHOLD
39 MERGE_THRESH_ZDIST = 30.0; // FIBER DISTANCE THRESHOLD
40 MERGE_THRESH_XYDIST = 40.0; // CENTROID TRANSVERSE DISTANCE THRESHOLD
41
42 // this parameter is used to break clusters based on rms time
43 BREAK_THRESH_TRMS= 5.0; // T RMS THRESHOLD
44
45 gPARMS->SetDefaultParameter( "BCALRECON:CLUST_THRESH", CLUST_THRESH );
46 gPARMS->SetDefaultParameter( "BCALRECON:MERGE_THRESH_DIST", MERGE_THRESH_DIST );
47 gPARMS->SetDefaultParameter( "BCALRECON:MERGE_THRESH_TIME", MERGE_THRESH_TIME );
48 gPARMS->SetDefaultParameter( "BCALRECON:MERGE_THRESH_ZDIST", MERGE_THRESH_ZDIST );
49 gPARMS->SetDefaultParameter( "BCALRECON:MERGE_THRESH_XYDIST", MERGE_THRESH_XYDIST );
50 gPARMS->SetDefaultParameter( "BCALRECON:BREAK_THRESH_TRMS", BREAK_THRESH_TRMS );
51
52 if( !DBCALGeometry::summingOn() ){
53
54 // these are energy calibration parameters -- no summing of cells
55
56 m_scaleZ_p0 = 0.9597;
57 m_scaleZ_p1 = 0.000454875;
58 m_scaleZ_p2 = -2.29912e-06;
59 m_scaleZ_p3 = 1.49757e-09;
60
61 m_nonlinZ_p0 = -0.00154122;
62 m_nonlinZ_p1 = 6.73594e-05;
63 m_nonlinZ_p2 = 0;
64 m_nonlinZ_p3 = 0;
65
66 }
67 else{
68
69 // these are energy calibration parameters -- 1.2.3.4 summing
70
71 //last updated for svn revision 9233
72 m_scaleZ_p0 = 0.99798;
73 m_scaleZ_p1 = 0.000361096;
74 m_scaleZ_p2 = -2.17338e-06;
75 m_scaleZ_p3 = 1.32201e-09;
76
77 m_nonlinZ_p0 = -0.0201272;
78 m_nonlinZ_p1 = 0.000103649;
79 m_nonlinZ_p2 = 0;
80 m_nonlinZ_p3 = 0;
81 }
82
83}
84
85//------------------
86// brun
87//------------------
88jerror_t DBCALShower_factory_KLOE::brun(JEventLoop *loop, int runnumber)
89{
90
91 vector<const DBCALGeometry*> bcalGeomVect;
92 loop->Get( bcalGeomVect );
93 const DBCALGeometry& bcalGeom = *(bcalGeomVect[0]);
94
95 //////////////////////////////////////////////////////////////////
96 // Calculate Cell Position
97 //////////////////////////////////////////////////////////////////
98
99 ATTEN_LENGTH = bcalGeom.ATTEN_LENGTH;
100 C_EFFECTIVE = bcalGeom.C_EFFECTIVE;
101
102 fiberLength = bcalGeom.BCALFIBERLENGTH; // fiber length in cm
103 zOffset = bcalGeom.GLOBAL_CENTER;
104
105 //the following uses some sad notation in which modmin=0 and modmax=48, when in fact there are 48 modules labelled either 0-47 or 1-48 depending on one's whim, although if we are using the methods from DBCALGeometry (e.g. cellId()), we must start counting from 1 and if we are accessing arrays we must of course start from 0
106 int modmin = 0;
107 int modmax = bcalGeom.NBCALMODS;
108 int rowmin1=0;
109 int rowmax1= bcalGeom.NBCALLAYSIN;
110 int rowmin2= rowmax1;
111 int rowmax2= bcalGeom.NBCALLAYSOUT+rowmin2;
112 int colmin1=0;
113 int colmax1=bcalGeom.NBCALSECSIN;
114 int colmin2=0;
115 int colmax2=bcalGeom.NBCALSECSOUT;
116
117 float r_inner= bcalGeom.BCALINNERRAD;
118
119 for (int i = (rowmin1+1); i < (rowmax1+1); i++){
120 //this loop starts from 1, so we can use i in cellId with no adjustment
121 int cellId = bcalGeom.cellId(1,i,1); //this gives us a cellId that we can use to get the radius. the module and sector numbers are irrelevant
122 //rt is radius of center of layer - BCAL inner radius
123 rt[i]=bcalGeom.r(cellId)-r_inner;
124 }
125
126 for (int i = (rowmin2+1); i < (rowmax2+1); i++){
127 int cellId = bcalGeom.cellId(1,i,1);
128 rt[i]=bcalGeom.r(cellId)-r_inner;
129 }
130
131 //these are r and phi positions of readout cells
132 float r[modulemax_bcal48][layermax_bcal10][colmax_bcal4];
133 float phi[modulemax_bcal48][layermax_bcal10][colmax_bcal4];
134
135 // Now start to extract cell position information from Geometry class
136 for (int k = modmin; k < modmax; k++){
137 for (int i = rowmin1; i < rowmax1; i++){
138 for (int j = colmin1; j < colmax1; j++){
139 //in this case the loops start at 0, so we have to add 1 to the indices when calling cellId(). hooray!
140 //use DBCALGeometry to get r/phi position of each cell
141 int cellId = bcalGeom.cellId(k+1,i+1,j+1);
142 r[k][i][j]=bcalGeom.r(cellId);
143 phi[k][i][j]=bcalGeom.phi(cellId);
144 //set x and y positions
145 xx[k][i][j]=r[k][i][j]*cos(phi[k][i][j]);
146 yy[k][i][j]=r[k][i][j]*sin(phi[k][i][j]);
147 }
148 }
149
150 for (int i = rowmin2; i < rowmax2; i++){
151 for (int j = colmin2; j < colmax2; j++){
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
162 ////////////////////////////////////////////////////////////////////////////
163 // Now the cell information are already contained in xx and yy arrays.
164 // xx and yy arrays are private members of this class
165 ////////////////////////////////////////////////////////////////////////////
166
167 // these are timing offsets -- a global side offset and individual cell offsets
168 ta0=0.0;
169 tb0=0.0;
170
171 for (int i = 0; i < modulemax_bcal48; i++){
172 for (int j = 0; j < layermax_bcal10; j++){
173 for (int k = 0; k < colmax_bcal4; k++){
174
175 ta_offset[i][j][k]=0.0;
176 tb_offset[i][j][k]=0.0;
177 }
178 }
179 }
180
181 return NOERROR;
182}
183
184//------------------
185// evnt
186//------------------
187jerror_t DBCALShower_factory_KLOE::evnt(JEventLoop *loop, int eventnumber)
188{
189 // Needed for associated objects later
190 vector<const DBCALHit*> bcalhits;
191 loop->Get(bcalhits);
192
193 // Call core KLOE reconstruction routines
194 CellRecon(loop);
195 CeleToArray();
196 PreCluster(loop);
197 ClusNorm();
198 ClusAnalysis();
199 Trakfit();
200
201 // Loop over reconstructed clusters and make DBCALShower objects out of them
202 vector<DBCALShower*> clusters;
203
204 int id = 0;
205 for (int i = 1; i < (clstot+1); i++){
206
207 int j=clspoi[i];
208
209 if( e_cls[j] < CLUST_THRESH ) continue;
210
211 // Time to cook a final shower
212 DBCALShower *shower = new DBCALShower;
213
214 shower->id = id++;
215 shower->E_raw = e_cls[j];
216 shower->x = x_cls[j];
217 shower->y = y_cls[j];
218 shower->z = z_cls[j] + zOffset;
219 shower->t = t_cls[j];
220 shower->N_cell = ncltot[j];
221
222 shower->xErr = eapx[1][j];
223 shower->yErr = eapx[2][j];
224 shower->zErr = eapx[3][j];
225
226 shower->tErr = 0.5 * sqrt( trms_a[j] * trms_a[j] +
227 trms_b[j] * trms_b[j] );
228
229 // calibrate energy:
230 // Energy calibration has a z dependence -- the
231 // calibration comes from fitting E_rec / E_gen to scale * E_gen^nonlin
232 // for slices of z. These fit parameters (scale and nonlin) are then plotted
233 // as a function of z and fit.
234
235 // center of target should be a geometry lookup
236 float zTarget = 65*k_cm;
237 float r = sqrt( shower->x * shower->x + shower->y * shower->y );
238
239 float zEntry = ( shower->z - zTarget ) * ( DBCALGeometry::BCALINNERRAD / r );
240
241 float scale = m_scaleZ_p0 + m_scaleZ_p1*zEntry +
242 m_scaleZ_p2*(zEntry*zEntry) + m_scaleZ_p3*(zEntry*zEntry*zEntry);
243 float nonlin = m_nonlinZ_p0 + m_nonlinZ_p1*zEntry +
244 m_nonlinZ_p2*(zEntry*zEntry) + m_nonlinZ_p3*(zEntry*zEntry*zEntry);
245
246 shower->E = pow( (shower->E_raw ) / scale, 1 / ( 1 + nonlin ) );
247
248 // Trace back to the DBCALHit objects used in this shower and
249 // add them as associated objects.
250 vector<const DBCALHit*> hitsInShower;
251 FindHitsInShower(j, bcalhits, hitsInShower);
252 for(unsigned int j=0; j<hitsInShower.size(); j++){
253 shower->AddAssociatedObject(hitsInShower[j]);
254 }
255
256 _data.push_back(shower);
257 }
258
259 return NOERROR;
260}
261
262
263//------------------
264// FindHitsInShower()
265//------------------
266void DBCALShower_factory_KLOE::FindHitsInShower(int indx, vector<const DBCALHit*> &bcalhits, vector<const DBCALHit*> &hitsInShower)
267{
268 /// This is called after the clusters have been completely formed. Our
269 /// job is simply to find the DBCALHit objects used to form a given cluster.
270 /// This is so the DBCALHit objects can be added to the DBCALShower object
271 /// as AssociatedObjects.
272
273 // The variable indx indexes the next[] array as the starting cell for the
274 // cluster. It also indexes the narr[][] array which holds the module, layer,
275 // sector(column) values for the hits.
276 //
277 // Here, we need to loop over next[] elements starting at next[indx] until
278 // we find the one pointing to element "indx" (i.e. the start of the list of
279 // cells in the cluster.) For each of these, we must find the LAST
280 // member in bcalhits to have the same module, layer, number indicating
281 // that is a DBCALHit to be added.
282
283 int start_indx = indx;
284 do{
285 int module = narr[1][indx];
286 int layer = narr[2][indx];
287 int sector = narr[3][indx];
288
289 // Loop over BCAL hits, trying to find this one
290 const DBCALHit *uphit=NULL__null;
291 const DBCALHit *downhit=NULL__null;
292 for(unsigned int i=0; i<bcalhits.size(); i++){
293 if(bcalhits[i]->module !=module)continue;
294 if(bcalhits[i]->layer !=layer)continue;
295 if(bcalhits[i]->sector !=sector)continue;
296
297 if(bcalhits[i]->end == DBCALGeometry::kUpstream) uphit = bcalhits[i];
298 if(bcalhits[i]->end == DBCALGeometry::kDownstream) downhit = bcalhits[i];
299 }
300
301 if(uphit)hitsInShower.push_back(uphit);
302 if(downhit)hitsInShower.push_back(downhit);
303
304 indx = next[indx];
305 }while(indx != start_indx);
306
307}
308
309
310//------------------
311// CellRecon()
312//------------------
313void DBCALShower_factory_KLOE::CellRecon(JEventLoop *loop)
314{
315 //======================================================================
316 // This code is used to reconstruct cell information cell by cell.
317 // This code is adapted from kloe code by Chuncheng Xu. June 29,2005
318
319 //**********************************************************************
320 // The main purpose of this function is extracting information
321 // from DBCALHit
322 // objects to form the arrays: ecel_a,tcel_a,ecel_b,tcel_b
323 // and xcel,ycel,zcel,tcel,ecel,tcell_anor,tcell_bnor;
324 // Among these arrays,
325 // ecel_a,tcel_a,ecel_b,tcel_b, xcel,ycel,zcel,tcel,ecel,tcell_anor and
326 // tcell_bnor are the input of function CeleToArray()
327 // These arrys are 3-D arrays and are rather logically are indexed by
328 // [module][layer][column]
329 //**********************************************************************
330
331 /////////////////////////////////////////////////////////////////////
332 // Now start to take take out DBCALHit to form our private
333 // data member ecel_a,tcel_a,ecel_b,tcel_b
334 /////////////////////////////////////////////////////////////////////
335
336 memset( ecel_a, 0, modulemax_bcal48 * layermax_bcal10 *
337 colmax_bcal4 * sizeof( float ) );
338 memset( tcel_a, 0, modulemax_bcal48 * layermax_bcal10 *
339 colmax_bcal4 * sizeof( float ) );
340 memset( ecel_b, 0, modulemax_bcal48 * layermax_bcal10 *
341 colmax_bcal4 * sizeof( float ) );
342 memset( tcel_b, 0, modulemax_bcal48 * layermax_bcal10 *
343 colmax_bcal4 * sizeof( float ) );
344
345 // extract the BCAL hits
346
347 vector<const DBCALHit*> hits;
348 loop->Get(hits);
349 if(hits.size() <= 0) return;
350
351
352 for (unsigned int i = 0; i < hits.size(); i++) {
353
354 const DBCALHit *hit = hits[i];
355 int module = hit->module;
356 int layer = hit->layer;
357 int sector = hit->sector;
358 int end = hit->end;
359 float E = hit->E;
360 float t = hit->t;
361
362 // cache hits indexed by module/layer/sector
363 // this allows end to end matching in next step
364
365 if(end==0) { // upstream
366 ecel_a[module-1][layer-1][sector-1] = E;
367 tcel_a[module-1][layer-1][sector-1] = t;
368 }
369 else if(end==1) { // downstream
370 ecel_b[module-1][layer-1][sector-1] = E;
371 tcel_b[module-1][layer-1][sector-1] = t;
372 }
373 else
374 {
375 cout<<"no such end, it is neither side A Nor B \n";
376 }
377 }
378
379 ////////////////////////////////////////////////////////////////////////
380 // Now data member ecel_a,tcel_a,ecel_b,tcel_b are readily filled.
381 ////////////////////////////////////////////////////////////////////////
382
383 ////////////////////////////////////////////////////////////////////////
384 // Now start to reconstruct cell information from data member ecel_a,tcel_a
385 // ecel_b,tcel_b. The reconstructed cell information are stored in
386 // data member arrays xcel,ycel,zcel,tcel,ecel,tcell_anor,tcell_bnor.
387 ////////////////////////////////////////////////////////////////////////
388
389 //
390 // ************ Loop over all CELE elements ********************
391 //
392 // K module,I layer, J sector
393
394 for (int k = 0; k < modulemax_bcal48; k++){
395 for (int i = 0; i < layermax_bcal10; i++){
396 for (int j = 0; j < colmax_bcal4; j++){
397
398 // get times
399 float ta = tcel_a[k][i][j];
400 float tb = tcel_b[k][i][j];
401
402 if( ( ta == 0 ) || ( tb == 0 ) ||
403 ( fabs( ta - tb ) > 80. ) ) {
404
405 xcel[k][i][j] = 0.0;
406 ycel[k][i][j] = 0.0;
407 zcel[k][i][j] = 0.0;
408 tcel[k][i][j] = 0.0;
409 ecel[k][i][j] = 0.0;
410
411 continue;
412 }
413
414 float ea = ecel_a[k][i][j];
415 float eb = ecel_b[k][i][j];
416
417 // algorithm requires both ends be hit
418 // improve energy resolution with single hits?
419 if( ( ea < ethr_cell ) ||
420 ( eb < ethr_cell ) ){
421
422 xcel[k][i][j] = 0.0;
423 ycel[k][i][j] = 0.0;
424 zcel[k][i][j] = 0.0;
425 tcel[k][i][j] = 0.0;
426 ecel[k][i][j] = 0.0;
427 continue;
428 }
429
430
431 float ta_cal = ta - ta0 - ta_offset[k][i][j];
432 float tb_cal = tb - tb0 - tb_offset[k][i][j];
433
434 float x = xx[k][i][j];
435 float y = yy[k][i][j];
436
437 float z1 = -0.5 * fiberLength; // cm
438 float z2 = 0.5 * fiberLength; // cm
439 float z0 = 0.5 * ( z1 + z2 );
440
441 float z = 0.5 * C_EFFECTIVE * ( ta_cal - tb_cal ) + z0;
442
443 if( z > ( 0.5 * fiberLength ) ){
444
445 z = 0.5*fiberLength;
446 }
447 else if( z < ( -0.5 * fiberLength ) ){
448
449 z = -0.5*fiberLength;
450 }
451
452 float t = 0.5 * ( ta_cal + tb_cal - fiberLength / C_EFFECTIVE );
453
454 float dpma = min( ( ta_cal - t ) * C_EFFECTIVE, fiberLength );
455 float dpmb = min( ( tb_cal - t ) * C_EFFECTIVE, fiberLength );
456 float atta = exp( -dpma / ATTEN_LENGTH );
457 float attb = exp( -dpmb / ATTEN_LENGTH );
458
459 // could compute energy weighted average instead of straight average
460 float e = ( ( ea / atta ) + ( eb / attb ) ) / 2;
461
462 xcel[k][i][j] = x;
463 ycel[k][i][j] = y;
464 zcel[k][i][j] = z;
465 tcel[k][i][j] = t;
466 ecel[k][i][j] = e;
467 tcell_anor[k][i][j] = ta_cal;
468 tcell_bnor[k][i][j] = tb_cal;
469 }
470 }
471 }
472}
473
474
475//------------------
476// CeleToArray()
477//------------------
478void DBCALShower_factory_KLOE::CeleToArray(void)
479{
480 // THis code is adpapted from kloe code by Chuncheng Xu on June29,2005
481 // The following part is taken from kaloe's clurec_lib.f's subroutine
482 // cele_to_arr.
483 //
484 // It's original purpose is to extract the data information from array RW
485 // and IW into array NARR(j,i), CELDATA(j,i), nclus(i),next(i),
486 // and e_cel(i), x_cel(i),y_cel(i), z_cel(i), t_cel(i),ta_cel(i)
487 // and tb_cel(i)
488 // In the above explanationn, i is the index for cell, and different j
489 // is for different component for such cell.
490 // for example, NARR(1,i), NARR(2,i),NARR(3,i) are for module number,
491 // layer number and sector number in halld bcal simulation
492
493
494 // Now we assume that the information is stored in arrays ECEL_A(k,i,j)
495 // , ECEL_B(k,i,j), TCEL_A(k,i,j), TCEL_B(k,i,j), XCEL(k,I,J),YCEL(k,I,J),
496 // ZCEL(K,I,J), TCEL(K,I,J), ECEL(K,I,J),TCELL_ANOR(K,I,J),TCELL_BNOR(K,I,J),
497 // which are passed into here by event.inc through common block.
498
499 // these values e_cel(i), x_cel(i),y_cel(i), z_cel(i), t_cel(i),ta_cel(i)
500 // and tb_cel(i)) are extremly useful in later's clusterization
501 // and they are passed by common block to precluster subroutine too. (through
502 // clurec_cal.inc)
503 //
504 // we hope this is a good "bridge" from raw data to precluster.
505 // This way we can keep good use of their strategy of clusterization
506 // as much as possible, although we know that Fortran has it's bad fame
507 // of not so easy to be reused.
508
509
510 // Essentially this function takes the cell information (e.g. E,x,y,z,t)
511 // from the 3D arrays filled in CellRecon() and puts it into 1D arrays.
512 // These 1D arrays are indexed by an identifier. The module/layer/sector
513 // associated with this identifier can be found using the narr[][] array,
514 // as described above.
515
516 celtot=0;
517
518 for (int k = 0; k < modulemax_bcal48; k++){
519 for (int i = 0; i < layermax_bcal10; i++){
520 for (int j = 0; j < colmax_bcal4; j++){
521
522 float ea = ecel_a[k][i][j];
523 float eb = ecel_b[k][i][j];
524 float ta = tcel_a[k][i][j];
525 float tb = tcel_b[k][i][j];
526
527 if( (min(ea,eb)>ethr_cell) & (fabs(ta-tb)<35.) & (ta!=0.) & (tb!=0.)) {
528 celtot=celtot+1;
529 } else {
530 continue;
531 }
532
533
534 if(celtot>cellmax_bcal48*10*4) {
535 break;
536 }
537
538 narr[1][celtot]=k+1; // these numbers will
539 narr[2][celtot]=i+1; // be used by preclusters
540 narr[3][celtot]=j+1; // which will start from index of 1
541 // rather than from 0.
542
543 // why 0.145? -- these variables are used as weights
544 celdata[1][celtot]=ea/0.145;
545 celdata[2][celtot]=eb/0.145;
546
547 nclus[celtot] = celtot;
548 next[celtot] = celtot;
549
550 e_cel[celtot] = ecel[k][i][j];
551 x_cel[celtot] = xcel[k][i][j];
552 y_cel[celtot] = ycel[k][i][j];
553 z_cel[celtot] = zcel[k][i][j];
554 t_cel[celtot] = tcel[k][i][j];
555
556 ta_cel[celtot]=tcell_anor[k][i][j];
557 tb_cel[celtot]=tcell_bnor[k][i][j];
558 }
559 }
560 }
561}
562
563
564
565//------------------
566// PreCluster()
567//------------------
568void DBCALShower_factory_KLOE::PreCluster(JEventLoop *loop)
569{
570 //what this function does: for each cell with a hit it finds the maximum
571 //energy neighbor and Connect()'s the two. Two cells are neighbors if they
572 //are within one column of each other and within one row. The situation is
573 //slightly more complicated for two cells on opposite sides of the boundary
574 //between inner cells and outer and is described in more detail below, but
575 //essentially works out the same way. The purpose of Connect() is described
576 //in that function itself.
577
578 int k=1; // NUMBER OF NEARBY ROWS &/OR TO LOOK FOR MAX E CELL
579
580 // extract the BCAL Geometry
581 vector<const DBCALGeometry*> bcalGeomVect;
582 loop->Get( bcalGeomVect );
583 const DBCALGeometry& bcalGeom = *(bcalGeomVect[0]);
584
585 // calculate cell position
586
587 int modmin = 0;
588 int modmax = bcalGeom.NBCALMODS;
589
590
591 //these values make sense as actually minima/maxima if it is implied that rowmin1=1,colmin1=1,colmin2=1
592 int rowmax1= bcalGeom.NBCALLAYSIN;
593 int rowmin2= rowmax1+1;
594 //int rowmax2= bcalGeom.NBCALLAYSOUT+rowmin2-1;
595 int colmax1=bcalGeom.NBCALSECSIN;
596 int colmax2=bcalGeom.NBCALSECSOUT;
597
598 float r_middle= bcalGeom.BCALMIDRAD;
599
600 //radial size of the outermost inner layer
601 float thick_inner=bcalGeom.rSize(bcalGeom.cellId(1,bcalGeom.NBCALLAYSIN,1));
602 //radial size of the innermost outer layer
603 float thick_outer=bcalGeom.rSize(bcalGeom.cellId(1,bcalGeom.NBCALLAYSIN+1,1));
604
605 // this is the radial distance between the center of the innermost outer layer and the outermost inner layer
606 float dis_in_out=bcalGeom.r(bcalGeom.cellId(1,bcalGeom.NBCALLAYSIN+1,1))-bcalGeom.r(bcalGeom.cellId(1,bcalGeom.NBCALLAYSIN,1));
607
608 float degree_permodule=360.0/(modmax-modmin);
609 float half_degree_permodule=degree_permodule/2.0;
610
611 //roughly the width of a single cell in the outermost inner layer
612 float width_1=2.0*(r_middle-thick_inner/2.0)*
613 sin(half_degree_permodule*3.141593/180)/colmax1;
614 //roughly the width of a single cell in the innermost outer layer
615 float width_2=2.0*(r_middle+thick_outer/2.0)*
616 sin(half_degree_permodule*3.141593/180)/colmax2;
617
618 //disthres is roughly the azimuthal distance between the center of an outer cell and the center of the most distant inner cell bordering an adjacent outer cell (a picture would be nice wouldn't it)
619 //this value is used for determining if two cells that straddle the boundary between inner and outer layers should be considered as neighboring
620 float disthres=width_2*1.5-width_1*0.5+0.0001;
621
622 for (int i = 1; i < (celtot+1); i++){
623
624 int maxnn=0; //cell index of the maximum energy neighbor (if one is found)
625 float emin=0.; //energy of maximum energy neighbor
626
627
628 for (int j = 1; j < (celtot+1); j++){
629 if ( (j!=i) & (nclus[j]!=nclus[i]) & (e_cel[j]>emin)) {
630
631 int k1= narr[1][i];
632 int k2= narr[1][j];
633 int i1= narr[2][i];
634 int i2= narr[2][j];
635
636
637 int modiff = k1-k2;
638 int amodif = abs(modiff);
639
640 // the following if is to check module and row distance.
641 if ( (abs(i1-i2)<=k) & ((amodif<=1) || (amodif==47)) ) {
642 // further check col distance
643 int j1= narr[3][i];
644 int j2= narr[3][j];
645
646 if(amodif==0) { // same module
647 // further check col distance if both are inner layers
648 if ( (i1<=rowmax1) & (i2<=rowmax1) & (abs(j2-j1)<=k) ) {
649 emin=e_cel[j];
650 maxnn=j;
651 }
652
653 // further check col distance if both are outer layers
654
655 if ( (i1>=rowmin2) & (i2>=rowmin2) & (abs(j2-j1)<=k) ) {
656 emin=e_cel[j];
657 maxnn=j;
658 }
659 }
660
661 if(amodif>0) { // different module
662 if( (modiff==1) || (modiff==-47) ) {
663 if ( (i1<=rowmax1) & (i2<=rowmax1) ){
664 if(abs((j1+colmax1)-j2)<=k){
665 emin=e_cel[j];
666 maxnn=j;
667 }
668 }
669
670 if ( (i1>=rowmin2) & (i2>=rowmin2) ) {
671 if(abs((j1+colmax2)-j2)<=k){
672 emin=e_cel[j];
673 maxnn=j;
674 }
675 }
676 }
677
678 if ( (modiff==-1) || (modiff==47) ) {
679
680 if ( (i1<=rowmax1) & (i2<=rowmax1) ){
681 if(abs((j2+colmax1)-j1)<=k){
682 emin=e_cel[j];
683 maxnn=j;
684 }
685 }
686
687 if ( (i1>=rowmin2) & (i2>=rowmin2) ){
688 if(abs((j2+colmax2)-j1)<=k){
689 emin=e_cel[j];
690 maxnn=j;
691 }
692 }
693 }
694 }
695
696 // further check col distance if one is inner layer, another is outer
697 // so that the two may be between the boundary of two different size
698 // of cells.
699 if( ( (i1 == rowmax1) & (i2 == rowmin2) ) ||
700 ( (i1 == rowmin2) & (i2 == rowmax1) ) ) {
701
702 float delta_xx=xx[k1-1][i1-1][j1-1]-xx[k2-1][i2-1][j2-1];
703 float delta_yy=yy[k1-1][i1-1][j1-1]-yy[k2-1][i2-1][j2-1];
704
705 //distance between centers of two cells
706 float dis = sqrt( delta_xx * delta_xx + delta_yy * delta_yy );
707 //dis_in_out is the distance in radial direction, so we now isolate distance in direction perpendicular to radius
708 dis = sqrt( dis*dis - dis_in_out * dis_in_out );
709 //disthres is described above
710 if( dis < disthres ){
711 emin = e_cel[j];
712 maxnn = j;
713 }
714 }
715 }
716 }
717 } // finish second loop
718
719 if(maxnn>0){
720
721 Connect(maxnn,i);
722 }
723 } // finish first loop
724}
725
726
727
728//------------------
729// Connect(n,m);
730//------------------
731void DBCALShower_factory_KLOE::Connect(int n,int m)
732{
733 //----------------------------------------------------------------------
734 // Purpose and Methods :CONNECTS CELL M TO THE NEAREST MAX NEIGHBOR N
735 // Created 31-JUL-1992 WON KIM
736 //----------------------------------------------------------------------
737
738 // This little piece of code wasn't so easy to decipher, but I *think* I
739 // understand what it's doing. The idea is the following:
740 //
741 // (Prior to entering this routine)
742 // The sparsified list of hits is copied into some 1-D arrays with
743 // dimension cellmax_bcal+1. This includes the nclus[] and next[]
744 // arrays. The nclus[] array keeps the cluster number which is initalized
745 // to the sparsified hit number. Thus, every (double-ended) hit cell
746 // is it's own cluster. A loop over pairs of hits is done above to find
747 // nearest neighbors that should be merged into the same cluster.
748 //
749 // The next[] array contains an index to the "next" cell in the cluster
750 // for the given hit cell. This is a circular list such that if one follows
751 // the "next[next[next[...]]]" values, the last element points back to
752 // the first element of the cluster. As such, these are also initialized
753 // to index themselves as all single hits are considered a cluster of
754 // one element prior to calling this Connect() routine.
755 //
756 // When we are called, the value of "m" will be less than than value
757 // of "n". The cluster number (kept in nclus[]) is therefore updated
758 // for all members of nclus[m] to be the same as nclus[n]. Furthermore,
759 // the hit cell "m" is appended to the cluster nclus[n]. This also means
760 // that the cluster numbers become non-sequential since all elements of
761 // nclus whose value is "m" will be changed such that their values are
762 // "n".
763 //
764 // 5/10/2011 DL
765
766 if(nclus[n]!=nclus[m]){
767 int j=m;
768 nclus[j]=nclus[n];
769 while(next[j]!=m){
770 j=next[j];
771 nclus[j]=nclus[n];
772 }
773 next[j]=next[n];
774 next[n]=m;
775 }
776}
777
778
779//------------------
780// ClusNorm()
781//------------------
782void DBCALShower_factory_KLOE::ClusNorm(void)
783{
784 // fast initialization of arrays:
785 memset( e_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
786 memset( x_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
787 memset( y_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
788 memset( z_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
789 memset( t_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
790 memset( ea_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
791 memset( eb_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
792 memset( ta_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
793 memset( tb_cls, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
794 memset( tsqr_a, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
795 memset( tsqr_b, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
796 memset( trms_a, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
797 memset( trms_b, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
798 memset( e2_a, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
799 memset( e2_b, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
800 memset( clspoi, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
801 memset( ncltot, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
802 memset( ntopol, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
803
804 // Part of what is being done here is to further sparsify the
805 // data into a list of clusters. This starts to fill arrays
806 // with dimension clsmax_bcal+1. One important thing is that
807 // the clspoi[] array is being filled with the cluster number
808 // of what should be the index of the first cell hit in the
809 // cluster.
810 //
811 // 5/10/2011 DL
812
813 clstot=0;
814
815 for (int ix = 1; ix < (celtot+1); ix++){
816
817 //----------------------------------------------------------------------
818 // KEEP TALLY OF CLUSTERS
819 //----------------------------------------------------------------------
820
821 int n=nclus[ix];
822 int j=0;
823
824 for (int i = 1; i < (clstot+1); i++){
825 if(n==clspoi[i]) j=i;
826 }
827
828 if(j==0) {
829 clstot=clstot+1;
830 clspoi[clstot]=n;
831 }
832
833 //----------------------------------------------------------------------
834 // NORMALIZED QUANTITIES OF THE TOTAL CLUSTER
835 //----------------------------------------------------------------------
836 if(e_cel[ix]<0.000000001)continue; // just for protection
837
838 x_cls[n]=(e_cls[n]*x_cls[n]+e_cel[ix]*x_cel[ix])
839 /(e_cls[n]+e_cel[ix]);
840
841 y_cls[n]=(e_cls[n]*y_cls[n]+e_cel[ix]*y_cel[ix])
842 /(e_cls[n]+e_cel[ix]);
843
844 z_cls[n]=(e_cls[n]*z_cls[n]+e_cel[ix]*z_cel[ix])
845 /(e_cls[n]+e_cel[ix]);
846
847 t_cls[n]=(e_cls[n]*t_cls[n]+e_cel[ix]*t_cel[ix])
848 /(e_cls[n]+e_cel[ix]);
849
850 e_cls[n]=e_cls[n]+e_cel[ix];
851
852 // write(*,*)' x-y-z-T-e done'
853 //----------------------------------------------------------------------
854 // NORMALIZED QUANTITIES FOR EACH SIDE
855 //----------------------------------------------------------------------
856
857 ta_cls[n]=(ea_cls[n]*ta_cls[n]+celdata[1][ix]*ta_cel[ix])
858 /(ea_cls[n]+celdata[1][ix]);
859 tsqr_a[n]=(ea_cls[n]*tsqr_a[n]+celdata[1][ix]*ta_cel[ix]*
860 ta_cel[ix])/(ea_cls[n]+celdata[1][ix]);
861 ea_cls[n]=ea_cls[n]+celdata[1][ix];
862 e2_a[n]=e2_a[n]+celdata[1][ix]*celdata[1][ix];
863 tb_cls[n]=(eb_cls[n]*tb_cls[n]+celdata[2][ix]*tb_cel[ix])/
864 (eb_cls[n]+celdata[2][ix]);
865 tsqr_b[n]=(eb_cls[n]*tsqr_b[n]+celdata[2][ix]*tb_cel[ix]*
866 tb_cel[ix])/(eb_cls[n]+celdata[2][ix]);
867 eb_cls[n]=eb_cls[n]+celdata[2][ix];
868 e2_b[n]=e2_b[n]+celdata[2][ix]*celdata[2][ix];
869
870 //----------------------------------------------------------------------
871 // TOTAL CELLS AND CELLS IN OTHER MODULES
872 //----------------------------------------------------------------------
873
874 ncltot[n]++;
875
876 if( narr[1][n] != narr[1][ix] || narr[2][n] != narr[2][ix] )
877 ntopol[n]++;
878
879 //----------------------------------------------------------------------
880
881 }
882
883 for (int n = 1; n < (clstot+1); n++){
884
885 int ix=clspoi[n];
886 if( ncltot[ix] > 1) {
887
888 float effnum = ea_cls[ix] * ea_cls[ix] / e2_a[ix];
889 trms_a[ix] = effnum / ( effnum - 1 ) *
890 ( tsqr_a[ix] - ta_cls[ix] * ta_cls[ix] );
891
892 effnum = eb_cls[ix] * eb_cls[ix] / e2_b[ix];
893 trms_b[ix] = effnum / ( effnum - 1 ) *
894 ( tsqr_b[ix] - tb_cls[ix] * tb_cls[ix] );
895
896 if( trms_a[ix] <= 0.0 ) trms_a[ix] = 0.;
897 if( trms_b[ix] <= 0.0 ) trms_b[ix] = 0.;
898 trms_a[ix] = sqrt( trms_a[ix] );
899 trms_b[ix] = sqrt( trms_b[ix] );
900 }
901 else {
902 trms_a[ix] = 0.;
903 trms_b[ix] = 0.;
904 }
905 }
906}
907
908//------------------
909// ClusAnalysis()
910//------------------
911void DBCALShower_factory_KLOE::ClusAnalysis()
912{
913 // track when clusters change to cut down
914 // on excess calles to expensive ClusNorm
915 bool newClust = false;
916
917 //----------------------------------------------------------------------
918 // Check for overlapping clusters
919 //----------------------------------------------------------------------
920
921 for (int i = 0; i < 2; i++){
922 for (int j = 1; j < (clstot+1); j++){
923 int ix=clspoi[j];
924 if(e_cls[ix]>0.0){
925 float dist=sqrt(trms_a[ix]*trms_a[ix]+trms_b[ix]*trms_b[ix]);
926 if(dist>BREAK_THRESH_TRMS) {
927 Clus_Break(ix);
928 newClust = true;
929 }
930 }
931 }
932
933 if( newClust ){
934
935 ClusNorm();
936 newClust = false;
937 }
938 }
939
940 //----------------------------------------------------------------------
941 // merge clusters likely to be from the same shower
942 //----------------------------------------------------------------------
943
944 int icls[3];
945 for (int i = 1; i < clstot; i++){
946 icls[1]=0;
947 icls[2]=0;
948 for (int j = (i+1); j < (clstot+1); j++){
949
950 int ix=clspoi[i];
951 int iy=clspoi[j];
952
953
954 if ( (e_cls[ix]>0.0) & (e_cls[iy]>0.0) ) {
955
956 float delta_x=x_cls[ix]-x_cls[iy];
957 float delta_y=y_cls[ix]-y_cls[iy];
958 float delta_z=z_cls[ix]-z_cls[iy];
959 float dist=sqrt(delta_x*delta_x+delta_y*delta_y+delta_z*delta_z);
960
961 float tdif=fabs(t_cls[ix]-t_cls[iy]);
962
963 // float distz=abs(z_cls[ix]-z_cls[iy]);
964 // cout<<"dist="<<dist<<" distz="<<distz<<" tdif="<<tdif<<"\n";
965
966 if ( (dist<MERGE_THRESH_DIST) & (tdif<MERGE_THRESH_TIME) ){
967 float zdif=fabs(z_cls[ix]-z_cls[iy]);
968 float distran=sqrt(delta_x*delta_x+delta_y*delta_y);
969
970 if ( (zdif<MERGE_THRESH_ZDIST) & (distran<MERGE_THRESH_XYDIST) ){
971 if(e_cls[ix]>=e_cls[iy]) {
972 icls[1]=ix;
973 icls[2]=iy;
974 }
975 else {
976 icls[1]=iy;
977 icls[2]=ix;
978 }
979 }
980 }
981
982 }
983
984 if(min(icls[1],icls[2])>0){
985
986 Connect(icls[1],icls[2]);
987 newClust = true;
988 }
989 }
990 }
991
992 if( newClust ){
993
994 ClusNorm();
995 }
996}
997
998//------------------
999// Clus_Break(ix);
1000//------------------
1001void DBCALShower_factory_KLOE::Clus_Break(int nclust)
1002{
1003 int nseed[5],selnum,selcel[cellmax_bcal48*10*4+1];
1004 float tdif,tdif_a,tdif_b,tseed[5];
1005
1006 //----------------------------------------------------------------------
1007 for (int i =0; i < 5; i++){
1008 nseed[i]=0;
1009 tseed[i]=0;
1010 }
1011 //----------------------------------------------------------------------
1012 // Divide cluster cells into four quadrant groups
1013 //----------------------------------------------------------------------
1014 int n=nclust;
1015 tdif_a=ta_cel[n]-ta_cls[nclust];
1016 tdif_b=tb_cel[n]-tb_cls[nclust];
1017 selnum=0;
1018
1019 //----------------------------------------------------------------------
1020 if(tdif_a>0.0) {
1021 if(tdif_b>0){
1022 selnum=1;
1023 }
1024 else {
1025 selnum=2;
1026 }
1027 }
1028
1029 else {
1030 if(tdif_b>0.0) {
1031 selnum=3;
1032 }
1033 else {
1034 selnum=4;
1035 }
1036 }
1037
1038 //------------------------------------------------------------------------
1039
1040 if(selnum>0) {
1041 float tdif=sqrt(tdif_a*tdif_a+tdif_b*tdif_b);
1042 if(tdif>tseed[selnum]){
1043 nseed[selnum]=n;
1044 tseed[selnum]=tdif;
1045 }
1046 selcel[n]=selnum;
1047 }
1048
1049 //-------------------------------------------------------------------------
1050
1051 while(next[n]!=nclust) {
1052 n=next[n];
1053 tdif_a=ta_cel[n]-ta_cls[nclust];
1054 tdif_b=tb_cel[n]-tb_cls[nclust];
1055 selnum=0;
1056
1057 //**************************************************************************
1058
1059 if(tdif_a>0.0) {
1060
1061 if(tdif_b>0.0) {
1062 selnum=1;
1063 }
1064 else {
1065 selnum=2;
1066 }
1067 }
1068
1069 else {
1070 if(tdif_b>0.0) {
1071 selnum=3;
1072 }
1073 else {
1074 selnum=4;
1075 }
1076 }
1077
1078 //***************************************************************************
1079
1080 //...........................................................................
1081
1082
1083 if(selnum>0){
1084 tdif=sqrt(tdif_a*tdif_a+tdif_b*tdif_b);
1085
1086 if(tdif>tseed[selnum]){
1087 nseed[selnum]=n;
1088 tseed[selnum]=tdif;
1089 }
1090
1091 selcel[n]=selnum;
1092 }
1093
1094 //...........................................................................
1095
1096 } //this bracket is related to "while" above.
1097
1098
1099 //----------------------------------------------------------------------
1100 // If successful, divide cluster chain into the new cluster chains
1101 //----------------------------------------------------------------------
1102
1103 for (int i =1; i < 5; i++){
1104
1105 if(nseed[i]>0) {
1106
1107 nclus[nseed[i]]=nseed[i];
1108 next[nseed[i]]=nseed[i];
1109
1110 for (int j =1; j < (celtot+1); j++){
1111 if ( (nclus[j]==nclust) & (j!=nseed[i]) ){
1112 if(selcel[j]==i) {
1113 nclus[j]=j;
1114 next[j]=j;
1115 Connect(nseed[i],j);
1116 }
1117 }
1118 }
1119 }
1120 }
1121}
1122
1123
1124//------------------
1125// Trakfit()
1126//------------------
1127
1128// routine modified by MRS to do expensive filling of layer information
1129// contains code previously in ClusNorm which is called multiple times
1130// per event, but there is no need to call Trakfit multiple times per event
1131void DBCALShower_factory_KLOE::Trakfit( void )
1132{
1133
1134 float emin=0.0001;
1135
1136 memset( clslyr, 0, ( clsmax_bcal48*10*4 + 1 ) *
1137 ( layermax_bcal10 + 1 ) * 6 * sizeof( float ) );
1138
1139 memset( apx, 0, ( clsmax_bcal48*10*4 + 1 ) * 4 * sizeof( float ) );
1140 memset( eapx, 0, ( clsmax_bcal48*10*4 + 1 ) * 4 * sizeof( float ) );
1141 memset( ctrk, 0, ( clsmax_bcal48*10*4 + 1 ) * 4 * sizeof( float ) );
1142 memset( ectrk, 0, ( clsmax_bcal48*10*4 + 1 ) * 4 * sizeof( float ) );
1143
1144 for (int ix = 1; ix < (celtot+1); ix++){
1145
1146 int n = nclus[ix];
1147
1148 //----------------------------------------------------------------------
1149 // NORMALIZED QUANTITIES OF THE TOTAL CLUSTER PER LAYER
1150 //----------------------------------------------------------------------
1151
1152 int lyr = narr[2][ix];
1153
1154 // write(*,*)'X, E',E_CEL(ix),CLSLYR(XLYR,LYR,N)
1155 clslyr[xlyr][lyr][n]=(clslyr[elyr][lyr][n]*clslyr[xlyr][lyr][n]
1156 +e_cel[ix]*x_cel[ix])/(clslyr[elyr][lyr][n]+e_cel[ix]);
1157
1158 clslyr[ylyr][lyr][n]=(clslyr[elyr][lyr][n]*clslyr[ylyr][lyr][n]
1159 +e_cel[ix]*y_cel[ix])/(clslyr[elyr][lyr][n]+e_cel[ix]);
1160 // write(*,*)' Y'
1161
1162 clslyr[zlyr][lyr][n]=(clslyr[elyr][lyr][n]*clslyr[zlyr][lyr][n]
1163 +e_cel[ix]*z_cel[ix])/(clslyr[elyr][lyr][n]+e_cel[ix]);
1164
1165 // write(*,*)' z'
1166 clslyr[tlyr][lyr][n]=(clslyr[elyr][lyr][n]*clslyr[tlyr][lyr][n]
1167 +e_cel[ix]*t_cel[ix])/(clslyr[elyr][lyr][n]+e_cel[ix]);
1168
1169 // write(*,*)'E'
1170 clslyr[elyr][lyr][n]=clslyr[elyr][lyr][n]+e_cel[ix];
1171
1172 // write(*,*)' slopes done'
1173 }
1174
1175 memset( nlrtot, 0, ( clsmax_bcal48*10*4 + 1 ) * sizeof( float ) );
1176
1177 for (int n = 1; n < ( clstot + 1 ); n++){
1178
1179 int ix=clspoi[n];
1180
1181 for (int i = 1; i < (layermax_bcal10+1); i++){
1182
1183 if( clslyr[elyr][i][ix] > 0.0 ) nlrtot[ix]++;
1184 }
1185
1186 for (int i = 0; i < ( layermax_bcal10 + 1 ); i++){
1187
1188 x[i]=0.0;
1189 y[i]=0.0;
1190 z[i]=0.0;
1191 e[i]=0.0;
1192 sigx[i]=0.0; // cm
1193 sigy[i]=0.0; // cm
1194 sigz[i]=0.0;
1195 }
1196
1197 int nltot=0;
1198
1199 for (int il = 1; il < (layermax_bcal10+1); il++){
1200
1201 if(clslyr[elyr][il][ix]>emin) {
1202
1203 nltot=nltot+1;
1204 x[nltot]= clslyr[xlyr][il][ix];
1205 y[nltot]= clslyr[ylyr][il][ix];
1206 z[nltot]= clslyr[zlyr][il][ix];
1207 e[nltot]= clslyr[elyr][il][ix];
1208
1209 sigy[nltot] = 1.0/e[nltot];
1210 sigx[nltot] = 1.0/e[nltot];
1211 sigz[nltot] = 1.0/sqrt(e[nltot]);
1212 }
1213 }
1214
1215 // The following error bar is the estimation of error bar
1216 // based on the experience of my fortran code running
1217 // If the structure of BCAL changes drastically, you have to make changes
1218 // accordingly. Xu Chuncheng , Jan 9, 2006
1219
1220 // **** this seems strange since it appears to overwrite what is above
1221
1222 sigx[1]=0.5;
1223 sigx[2]=0.5;
1224 sigx[3]=0.5;
1225 sigx[4]=0.5;
1226 sigx[5]=0.5;
1227 sigx[6]=0.8;
1228 sigx[7]=0.9;
1229 sigx[8]=1.2;
1230 sigx[9]=1.3;
1231
1232 sigy[1]=0.5;
1233 sigy[2]=0.5;
1234 sigy[3]=0.5;
1235 sigy[4]=0.5;
1236 sigy[5]=0.5;
1237 sigy[6]=0.8;
1238 sigy[7]=0.9;
1239 sigy[8]=1.2;
1240 sigy[9]=1.3;
1241
1242
1243 sigz[1]=0.5;
1244 sigz[2]=0.5;
1245 sigz[3]=0.5;
1246 sigz[4]=0.5;
1247 sigz[5]=0.5;
1248 sigz[6]=0.8;
1249 sigz[7]=0.9;
1250 sigz[8]=1.2;
1251 sigz[9]=1.3;
1252
1253 if( nltot > 1 ){
1254
1255 Fit_ls();
1256
1257 for (int i = 1; i < 4; i++){
1258
1259 ctrk[i][ix]=ctrk_ix[i];
1260 ectrk[i][ix]=ectrk_ix[i];
1261 apx[i][ix]=apx_ix[i];
1262 eapx[i][ix]=eapx_ix[i];
1263 }
1264 }
1265 else{
1266
1267 // we have 1 or less layers hit -- no fit
1268
1269 apx[1][ix] = x[1];
1270 apx[2][ix] = y[1];
1271 apx[3][ix] = z[1];
1272 eapx[1][ix] = sigx[1];
1273 eapx[2][ix] = sigy[1];
1274 eapx[3][ix] = sigz[1];
1275 ectrk[1][ix] = 0.0;
1276 ectrk[2][ix] = 0.0;
1277 ectrk[3][ix] = 0.0;
1278 ctrk[1][ix] = 999.0;
1279 ctrk[2][ix] = 999.0;
1280 ctrk[3][ix] = 999.0;
1281 }
1282 }
1283}
1284
1285//------------------
1286// Fit_ls()
1287//------------------
1288void DBCALShower_factory_KLOE::Fit_ls()
1289{
1290 float a,b,c;
1291 float d,e,f,chi2,q,norm;
1292 float siga,sigb,sigc,sigd,sige,sigf;
1293 float sigb2,sigd2,sigf2;
1294
1295 // fitting for X=a+bt
1296 Linefit(1,1,a,b,siga,sigb,chi2,q);
1297 // fitting for Y=c+dt
1298 Linefit(2,1,c,d,sigc,sigd,chi2,q);
1299 // fitting for Z=e+ft
1300 Linefit(3,1,e,f,sige,sigf,chi2,q);
1301 sigb2=sigb*sigb;
1302 sigd2=sigd*sigd;
1303 sigf2=sigf*sigf;
1304
1305 apx_ix[1]=a;
1306 apx_ix[2]=c;
1307 apx_ix[3]=e;
1308 eapx_ix[1]=siga;
1309 eapx_ix[2]=sigc;
1310 eapx_ix[3]=sige;
1311
1312 // write(*,*) "b,d,f = ",b,d,f
1313 // cout<<"b= "<<b<<" d="<<d<<" f="<<f<<"\n";
1314
1315 // write(*,*) "sigb,sigd,sigf = ",sigb,sigd,sigf
1316
1317 norm=sqrt(b*b+d*d+f*f);
1318
1319 ctrk_ix[1]=b/norm;
1320 ctrk_ix[2]=d/norm;
1321 ctrk_ix[3]=f/norm;
1322
1323 float norm3=norm*norm*norm;
1324
1325 ectrk_ix[1]=sqrt((d*d+f*f)*(d*d+f*f)*sigb2+b*b*d*d*sigd2+b*b*f*f*sigf2)/norm3;
1326 ectrk_ix[2]=sqrt((b*b+f*f)*(b*b+f*f)*sigd2+d*d*b*b*sigb2+d*d*f*f*sigf2)/norm3;
1327 ectrk_ix[3]=sqrt((b*b+d*d)*(b*b+d*d)*sigf2+f*f*b*b*sigb2+f*f*d*d*sigd2)/norm3;
1328
1329 return;
1330}
1331
1332
1333//------------------
1334// Linefit(x,y,ndata,sig,mwt,a,b,siga,sigb,chi2,q)
1335//------------------
1336void DBCALShower_factory_KLOE::Linefit(int ixyz,int mwt,float &a,
1337 float &b,float &siga,float &sigb,float &chi2,float &q)
1338{
1339
1340 // This programme is taken from Garth's book
1341 // "Numerical Recipes in Fortran"
1342 // and we tested it. It works well. Chuncheng Xu,2005
1343
1344
1345 float sig[layermax_bcal10+1],etemp;
1346 float xtemp[layermax_bcal10+1],ytemp[layermax_bcal10+1];
1347 // uses gammq
1348
1349 // Given a set of data points X(1:ndata),Y(1:ndata) with individual standard
1350 // deviations sig(1:ndata), fit them to a strait line y=a+bx by minimum
1351 // chisq. Returned are a,b and their respective probable uncertainties
1352 // siga and sigb, the chisq, and the goodness-of-fit probability q(that the
1353 // fit would have chisq this large or larger). if mwt=0 on input, then the
1354 // the standard deviations are assumed to be unavalable:q is returned as
1355 // 1.0 and the normalization of chi2 is to the unit standard deviation on all
1356 // points.
1357
1358
1359 float sigdat,ss,st2,sx,sxoss,sy,t,wt;
1360 sx=0.0; // Initialize sums to zero
1361 sy=0.0;
1362 st2=0.0;
1363 b=0.0;
1364
1365 int ndata=0;
1366
1367
1368
1369 if(ixyz==1) {
1370 for (int i = 1; i < (layermax_bcal10+1); i++){
1371 xtemp[i]=rt[i];
1372 ytemp[i]=x[i];
1373 sig[i]=sigx[i];
1374 etemp=e[i];
1375 if(etemp>0.0001)ndata=ndata+1;
1376
1377 }
1378 }
1379 else if(ixyz==2) {
1380 for (int i = 1; i < (layermax_bcal10+1); i++){
1381 xtemp[i]=rt[i];
1382 ytemp[i]=y[i];
1383 sig[i]=sigy[i];
1384 etemp=e[i];
1385 if(etemp>0.000001)ndata=ndata+1;
1386 }
1387 }
1388 else if(ixyz==3) {
1389 for (unsigned int i = 1; i < (layermax_bcal10+1); i++){
1390 xtemp[i]=rt[i];
1391 ytemp[i]=z[i];
1392 sig[i]=sigz[i];
1393 etemp=e[i];
1394 if(etemp>0.000001)ndata=ndata+1;
1395 }
1396 }
1397
1398 if(mwt!=0) { // Accumulate sums
1399 ss=0.0;
1400 for (int i = 1; i < (ndata+1); i++){
1401 wt=1.0/(sig[i]*sig[i]);
1402 ss=ss+wt;
1403 sx=sx+xtemp[i]*wt;
1404 sy=sy+ytemp[i]*wt;
1405 }
1406 }
1407
1408 else {
1409 for (int i = 1; i < (ndata+1); i++){
1410 sx=sx+xtemp[i];
1411 sy=sy+ytemp[i];
1412 }
1413 ss=float(ndata);
1414 }
1415
1416 sxoss=sx/ss;
1417
1418 if(mwt!=0) {
1419 for (int i = 1; i < (ndata+1); i++){
1420 t=(xtemp[i]-sxoss)/sig[i];
1421 st2=st2+t*t;
1422 b=b+t*ytemp[i]/sig[i];
1423 }
1424
1425 }
1426
1427 else {
1428
1429 for (int i = 1; i < (ndata+1); i++){
1430 t=xtemp[i]-sxoss;
1431 st2=st2+t*t;
1432 b=b+t*ytemp[i];
1433 }
1434 }
1435
1436 b=b/st2; // Solve for a,b,siga and sigb
1437 a=(sy-sx*b)/ss;
1438 siga=sqrt((1.0+sx*sx/(ss*st2))/ss);
1439 sigb=sqrt(1.0/st2);
1440 chi2=0.0; // Calculate chisq
1441
1442 if(mwt==0) {
1443 for (int i = 1; i < (ndata+1); i++){
1444 chi2=chi2+(ytemp[i]-a-b*xtemp[i])*(ytemp[i]-a-b*xtemp[i]);
1445 }
1446 q=1.0;
1447 sigdat=sqrt(chi2/(ndata-2));
1448 siga=siga*sigdat;
1449 sigb=sigb*sigdat;
1450 }
1451 else {
1452 for (int i = 1; i < (ndata+1); i++){
1453 chi2=chi2+((ytemp[i]-a-b*xtemp[i])/
1454 sig[i])*((ytemp[i]-a-b*xtemp[i])/sig[i]);
1455 }
1456 q=Gammq(0.5*(ndata-2),0.5*chi2);
1457 }
1458
1459}
1460
1461
1462
1463//------------------
1464// Gammq(a_gammq,x_gammq)
1465//------------------
1466float DBCALShower_factory_KLOE::Gammq(float a_gammq,float x_gammq)
1467{
1468
1469 //============================================================================
1470
1471
1472 float gammq;
1473
1474 // uses gcf,gser
1475 //Returns the incomplete gamma function Q(a_gammq,x_gammq)=1-P(a_gammq,x_gammq)
1476
1477
1478 float gammcf,gamser;
1
Variable 'gammcf' declared without an initial value
→
1479
1480 if(a_gammq==0.0) {
2
←
Taking false branch
→
1481 gammq=999.0;
1482 return gammq;
1483 }
1484
1485
1486
1487 if(x_gammq<0. || a_gammq<= 0.0) {
3
←
Taking false branch
→
1488 // cout<<"bad arguments in gammq"<<"\n";
1489 return 999.0;// pause 'bad arguments in gammq'
1490 }
1491
1492 if(x_gammq<(a_gammq+1.)) { // Use the series representation
4
←
Taking false branch
→
1493 Gser(gamser,a_gammq,x_gammq);
1494 gammq=1.0-gamser; // and take its complement
1495 }
1496 else { // Use continued fraction representation
1497 Gcf(gammcf,a_gammq,x_gammq);
5
←
Calling 'Gcf'
→
10
←
Returning from 'Gcf'
→
1498 gammq=gammcf;
11
←
Assigned value is garbage or undefined
1499 }
1500 return gammq;
1501}
1502
1503
1504//------------------
1505// Gser(gamser,a_gser,x_gser,gln)
1506//------------------
1507void DBCALShower_factory_KLOE::Gser(float &gamser,float a_gser,float x_gser)
1508{
1509 //============================================================================
1510 int itmax=100;
1511 float eps=3.0e-7;
1512 float gln;
1513
1514 // uses gammln
1515 // Returns the incomplete gamma function P(a,x) evaluated by its
1516 // series representation as gamser. Also returns ln(Gamma(a)) as gln
1517
1518 float ap, del,sum;
1519
1520 gln=Gammln(a_gser);
1521
1522 if(x_gser<=0.0) {
1523 if(x_gser<0.0) cout<<"x_gser<0 in gser"<<"\n";
1524 gamser=0.0;
1525 return;
1526 }
1527
1528 ap=a_gser;
1529 sum=1.0/a_gser;
1530 del=sum;
1531
1532
1533 for (int n = 1; n < (itmax+1); n++){
1534 ap=ap+1.0;
1535 del=del*x_gser/ap;
1536 sum=sum+del;
1537
1538 if(fabs(del)<fabs(sum)*eps) {
1539 gamser=sum*exp(-x_gser+a_gser*log(x_gser)-gln);
1540 return;
1541 }
1542
1543 }
1544
1545 // cout<< "a too large, ITMAX too small in gser"<<"\n";
1546 return;
1547}
1548
1549
1550//------------------
1551// Gcf(gammcf,a,x,gln)
1552//------------------
1553void DBCALShower_factory_KLOE::Gcf(float &gammcf,float a_gcf,float x_gcf)
1554{
1555 //========================================================================
1556
1557 int itmax=100;
1558 float eps=3.0e-7;
1559 float fpmin=1.0e-30;
1560
1561 float gln;
1562
1563 // uses gammln
1564 // Returns the incomplete gamma function Q(a,x) evaluated by
1565 // its continued fraction representation as gammcf. Also returns
1566 // ln(Gamma(a)) as gln
1567
1568 // Parameters: ITMAX is the maximum allowed number of iterations;
1569 // EPS is the relative accuracy; FPMIN is a number near the smallest
1570 // representable floating-point number.
1571
1572 float an,b,c,d,del,h;
1573
1574 gln=Gammln(a_gcf);
1575 b=x_gcf+1.0-a_gcf;
1576 c=1.0/fpmin;
1577 d=1.0/b;
1578 h=d;
1579
1580
1581 for (int i = 1; i < (itmax+1); i++){
6
←
Loop condition is true. Entering loop body
→
1582 an=-i*(i-a_gcf);
1583 b=b+2.0;
1584 d=an*d+b;
1585 if(fabs(d)<fpmin)d=fpmin;
7
←
Taking false branch
→
1586 c=b+an/c;
1587 if(fabs(c)<fpmin)c=fpmin;
8
←
Taking false branch
→
1588 d=1.0/d;
1589 del=d*c;
1590 h=h*del;
1591 if(fabs(del-1.0)<eps) {
9
←
Taking false branch
→
1592 gammcf=exp(-x_gcf+a_gcf*log(x_gcf)-gln)*h; // Put factors in front
1593 return;
1594 }
1595 // cout<< "a too large,ITMAX too small in gcf"<<"\n";
1596 return;
1597 }
1598} //???
1599
1600//------------------
1601// Gammln(xx)
1602//------------------
1603float DBCALShower_factory_KLOE::Gammln(float xx_gln)
1604{
1605 // Returns the value ln[Gamma(xx_gln)] for xx_gln>0.0
1606
1607 float ser,stp,tmp,x_gln,y_gln;
1608 float cof[7];
1609 float gammln;
1610
1611 // Internal arithmetic will be done in double precision, a nicety that
1612 // that you can omit if five-figure accuracy is good enoug
1613 // save cof,stp
1614 // data cof,stp/76.18009172947146d0,-86.50532032941677d0,
1615 // * 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
1616 // * -.5395239384953d-5,2.5066282746310005d0/
1617
1618 stp=2.5066282746310005;
1619 cof[1]=76.18009172947146;
1620 cof[2]=-86.50532032941677;
1621 cof[3]=24.01409824083091;
1622 cof[4]=-1.231739572450155;
1623 cof[5]=.1208650973866179e-2;
1624 cof[6]=-.5395239384953e-5;
1625
1626 x_gln=xx_gln;
1627 y_gln=x_gln;
1628 tmp=x_gln+5.5;
1629 tmp=(x_gln+0.5)*log(tmp)-tmp;
1630
1631 ser=1.000000000190015;
1632
1633 for (int j = 1; j < 7; j++){
1634 y_gln=y_gln+1.0;
1635 ser=ser+cof[j]/y_gln;
1636 }
1637
1638
1639 gammln=tmp+log(stp*ser/x_gln);
1640 return gammln;
1641}
1642