Bug Summary

File:/home/sdobbs/work/clang/halld_recon/src/libraries/BCAL/DBCALShower_factory_KLOE.cc
Warning:line 285, column 9
Value stored to 'sig_t' is never read

Annotated Source Code

Press '?' to see keyboard shortcuts

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