Bug Summary

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

Annotated Source Code

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