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 |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
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 | |
22 | using namespace std; |
23 | |
24 | //------------------ |
25 | // init |
26 | //------------------ |
27 | jerror_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 | //------------------ |
95 | jerror_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 | //------------------ |
184 | jerror_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 | //------------------ |
334 | void 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 | //------------------ |
379 | void 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 | //------------------ |
472 | void 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 | //------------------ |
562 | void 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 | //------------------ |
725 | void 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 | //------------------ |
776 | void 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 | //------------------ |
905 | void 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 | //------------------ |
995 | void DBCALShower_factory_KLOE::Clus_Break(int nclust) |
996 | { |
997 | int nseed[5],selnum,selcel[cellmax_bcal48*10*4+1]; |
998 | float tdif,tdif_a,tdif_b,tseed[5]; |
999 | |
1000 | //---------------------------------------------------------------------- |
1001 | for (int i =0; i < 5; i++){ |
1002 | nseed[i]=0; |
1003 | tseed[i]=0; |
1004 | } |
1005 | //---------------------------------------------------------------------- |
1006 | // 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 |
1125 | void DBCALShower_factory_KLOE::Trakfit( void ) |
1126 | { |
1127 | |
1128 | float emin=0.0001; |
1129 | |
1130 | memset( clslyr, 0, ( clsmax_bcal48*10*4 + 1 ) * |
1131 | ( layermax_bcal10 + 1 ) * 6 * sizeof( float ) ); |
1132 | |
1133 | memset( apx, 0, ( clsmax_bcal48*10*4 + 1 ) * 4 * sizeof( float ) ); |
1134 | memset( eapx, 0, ( clsmax_bcal48*10*4 + 1 ) * 4 * sizeof( float ) ); |
1135 | memset( ctrk, 0, ( clsmax_bcal48*10*4 + 1 ) * 4 * sizeof( float ) ); |
1136 | memset( ectrk, 0, ( clsmax_bcal48*10*4 + 1 ) * 4 * sizeof( float ) ); |
1137 | |
1138 | for (int ix = 1; ix < (celtot+1); ix++){ |
1139 | |
1140 | int n = nclus[ix]; |
1141 | |
1142 | //---------------------------------------------------------------------- |
1143 | // 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 | //------------------ |
1282 | void DBCALShower_factory_KLOE::Fit_ls() |
1283 | { |
1284 | float a,b,c; |
1285 | float d,e,f,chi2,q,norm; |
1286 | float siga,sigb,sigc,sigd,sige,sigf; |
1287 | float sigb2,sigd2,sigf2; |
1288 | |
1289 | // 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 | //------------------ |
1330 | void DBCALShower_factory_KLOE::Linefit(int ixyz,int mwt,float &a, |
1331 | float &b,float &siga,float &sigb,float &chi2,float &q) |
1332 | { |
1333 | |
1334 | // 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 | //------------------ |
1460 | float 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 | //------------------ |
1501 | void DBCALShower_factory_KLOE::Gser(float &gamser,float a_gser,float x_gser) |
1502 | { |
1503 | //============================================================================ |
1504 | int itmax=100; |
1505 | float eps=3.0e-7; |
1506 | float gln; |
1507 | |
1508 | // 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 | //------------------ |
1547 | void DBCALShower_factory_KLOE::Gcf(float &gammcf,float a_gcf,float x_gcf) |
1548 | { |
1549 | //======================================================================== |
1550 | |
1551 | int itmax=100; |
1552 | float eps=3.0e-7; |
1553 | float fpmin=1.0e-30; |
1554 | |
1555 | float gln; |
1556 | |
1557 | // 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 | //------------------ |
1597 | float 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 |