Bug Summary

File:libraries/BCAL/DBCALShower_factory_CURVATURE.cc
Location:line 376, column 5
Description:Value stored to 'sig_y' is never read

Annotated Source Code

1// $Id$
2//
3// File: DBCALShower_factory_CURVATURE.cc
4// Created: Fri Mar 27 10:57:45 CST 2015
5// Creator: beattite (on Linux eos.phys.uregina.ca 2.6.32-504.12.2.el6.x86_64 x86_64)
6//
7
8// This factory takes showers from the IU shower factory and alters them based on shower
9// curvature tables produced by AS. These can be found in the CCDB under
10// BCAL/curvature_central and BCAL/curvature_side. They give central z-positions of
11// energy depositions in the BCAL for each layer based on the reconstructed IU shower's
12// energy and theta values. The two tables are for 'central' sectors and 'side' sectors.
13// A central sector (the one or two sectors closest to the centroid of the shower) exhibit
14// different deposition ranges in z than a side sector.
15
16// We first grab the IU showers and check for overlap in z. We want to merge showers that
17// are close in z, t, and phi that the IU code may have reconstructed as two separate showers.
18// Next, we extract the energies and theta values for each of the IU (and merged) showers.
19// Using the curvature tables, we find the appropriate z-bin position for each layer, then
20// loop through the points in the event and put any point that falls within the z-bins for a
21// shower into that shower. The z-bin width is determined by the error on the z-positions of
22// the depositions and on a set width factor.
23
24// These new colletions of points are then averaged (energy-squared-weighted averages in x, y,
25// z, and t) to get the shower values.
26
27#include <iostream>
28#include <iomanip>
29using namespace std;
30
31#include "DBCALShower_factory_CURVATURE.h"
32#include "DBCALPoint.h"
33#include "DBCALShower_factory_IU.h"
34
35#include "DANA/DApplication.h"
36
37#include "units.h"
38using namespace jana;
39#include "TMath.h"
40
41//------------------
42// init
43//------------------
44jerror_t DBCALShower_factory_CURVATURE::init(void)
45{
46 if( ! DBCALGeometry::summingOn() ) {
47 // in libraries/PID/DNeutralShowerCandidate.h, there are some constants used to calculate the energy uncertainty. If you are updating these constants, you might want to update that also...
48
49 // these are energy calibration parameters -- no summing of cells
50
51 m_scaleZ_p0 = 0.950774;
52 m_scaleZ_p1 = 0.000483979;
53 m_scaleZ_p2 = -2.08086e-06;
54 m_scaleZ_p3 = 8.08534e-10;
55
56 m_nonlinZ_p0 = 0.0152548;
57 m_nonlinZ_p1 = 0;
58 m_nonlinZ_p2 = 0;
59 m_nonlinZ_p3 = 0;
60 }
61 else{
62
63 // these are energy calibration parameters -- 1.2.3.4 summing
64
65 //last updated for svn revision 9233
66 m_scaleZ_p0 = 0.992437;
67 m_scaleZ_p1 = 0.00039242;
68 m_scaleZ_p2 = -2.23135e-06;
69 m_scaleZ_p3 = 1.40158e-09;
70
71 m_nonlinZ_p0 = -0.0147086;
72 m_nonlinZ_p1 = 9.69207e-05;
73 m_nonlinZ_p2 = 0;
74 m_nonlinZ_p3 = 0;
75
76 }
77 return NOERROR;
78}
79
80//------------------
81// brun
82//------------------
83jerror_t DBCALShower_factory_CURVATURE::brun(jana::JEventLoop *loop, int32_t runnumber)
84{
85 DApplication* app = dynamic_cast<DApplication*>(loop->GetJApplication());
86 DGeometry* geom = app->GetDGeometry(runnumber);
87 geom->GetTargetZ(m_zTarget);
88
89 vector< vector<double> > curvature_parameters;
90
91 // Read in curvature parameters from two tables in the database.
92 // We'll use two four-dimentional arrays (position and sigma).
93 // The indices are: [sector type (central or side)][layer (from 0 to 4)][angle bin (from 0 to 11)][energy bin (from 0 to 31)]
94 for (int ii = 0; ii < 2; ii++){
95 if (ii == 0){
96 if (loop->GetCalib("/BCAL/curvature_central", curvature_parameters)) jout << "Error loading /BCAL/curvature_central !" << endl;
97 }
98 else {
99 if (loop->GetCalib("/BCAL/curvature_side", curvature_parameters)) jout << "Error loading /BCAL/curvature_side !" << endl;
100 }
101 for (int line = 0; line < 1536; line++){
102 layer = curvature_parameters[line][0];
103 angle = curvature_parameters[line][1];
104 energy = curvature_parameters[line][2];
105 dataposition = curvature_parameters[line][3];
106 datasigma = curvature_parameters[line][4];
107 if (angle == 115){ // angle bins are 5, 10, 20, 30, ..., 100, 110, 115. [angle/10 - 1] won't work for the 115 bin, so we explicitly set these.
108 position[ii][layer - 1][11][energy/50 - 1] = dataposition - 48; //table values measure from the end of the BCAL. Points measure from the center of the target.
109 sigma[ii][layer - 1][11][energy/50 - 1] = datasigma;
110 }
111 else {
112 position[ii][layer - 1][angle/10 - 1][energy/50 - 1] = dataposition - 48; //table values measure from the end of the BCAL. Points measure from the center of the target.
113 sigma[ii][layer - 1][angle/10 - 1][energy/50 - 1] = datasigma;
114 }
115 }
116 curvature_parameters.clear();
117 }
118
119 // Thresholds for merging showers and including points.
120 PHITHRESHOLD = 4*0.06544792; // 4*2 column widths, in radians.
121 ZTHRESHOLD = 35.*k_cm; // Range in z, in cm.
122 TTHRESHOLD = 8.0*k_nsec; // Loose time range, in ns.
123 ETHRESHOLD = 1.0*k_keV; // Points minimal energy, in keV.
124
125 return NOERROR;
126}
127
128//------------------
129// evnt
130//------------------
131jerror_t DBCALShower_factory_CURVATURE::evnt(JEventLoop *loop, uint64_t eventnumber)
132{
133 recon_showers_phi.clear();
134 recon_showers_theta.clear();
135 recon_showers_E.clear();
136 recon_showers_t.clear();
137 vector<const DBCALShower*> bcal_showers;
138 vector<const DBCALPoint*> Points;
139 vector<const DBCALPoint*> CurvaturePoints;
140
141 loop->Get(bcal_showers,"IU");
142 loop->Get(Points);
143
144// Merge overlapping showers and calculate new shower parameters (phi, theta, E)
145 while (bcal_showers.size() != 0){
146 double recon_x = bcal_showers[0]->x;
147 double recon_y = bcal_showers[0]->y;
148 double recon_phi = atan2(recon_y,recon_x);
149 double recon_theta = atan2(sqrt(bcal_showers[0]->x*bcal_showers[0]->x + bcal_showers[0]->y*bcal_showers[0]->y),bcal_showers[0]->z-m_zTarget);
150 double recon_E = bcal_showers[0]->E;
151 double recon_t = bcal_showers[0]->t/* - mcgentime*/;
152 double total_E = bcal_showers[0]->E;
153 for (int ii = 1; ii < (int)bcal_showers.size(); ii++){ // compare each shower in bcal_showers to the first to check for overlap.
154 double delPhi = atan2(bcal_showers[0]->y,bcal_showers[0]->x) - atan2(bcal_showers[ii]->y,bcal_showers[ii]->x);
155 if (delPhi > TMath::Pi()) delPhi -= 2*TMath::Pi();
156 if (delPhi < -1*TMath::Pi()) delPhi += 2*TMath::Pi();
157 double delZ = bcal_showers[0]->z - bcal_showers[ii]->z;
158 double delT = bcal_showers[0]->t - bcal_showers[ii]->t;
159 if (fabs(delPhi) < PHITHRESHOLD && fabs(delZ) < ZTHRESHOLD && fabs(delT) < TTHRESHOLD) overlap.push_back(ii); // Showers overlap!
160 }
161 if (overlap.size() != 0){
162 recon_x *= recon_E;
163 recon_y *= recon_E;
164 recon_theta *= recon_E;
165 recon_E *= recon_E;
166 recon_t *= recon_E;
167 }
168 while (overlap.size() != 0){ // If we had an overlap, average the overlapping showers together.
169 recon_x += bcal_showers[overlap.back()]->x*bcal_showers[overlap.back()]->E;
170 recon_y += bcal_showers[overlap.back()]->y*bcal_showers[overlap.back()]->E; // Average over x and y rather than over phi to avoid boundary errors.
171 recon_theta += atan2(sqrt(bcal_showers[overlap.back()]->x*bcal_showers[overlap.back()]->x + bcal_showers[overlap.back()]->y*bcal_showers[overlap.back()]->y),bcal_showers[overlap.back()]->z-m_zTarget)*bcal_showers[overlap.back()]->E;
172 recon_E += bcal_showers[overlap.back()]->E*bcal_showers[overlap.back()]->E;
173 recon_t += (bcal_showers[overlap.back()]->t/* - mcgentime*/)*bcal_showers[overlap.back()]->E;
174 total_E += bcal_showers[overlap.back()]->E;
175
176 bcal_showers.erase(bcal_showers.begin()+overlap.back()); // Erase the extra overlapping showers.
177 overlap.pop_back();
178 if (overlap.size() == 0){
179 recon_x /= total_E;
180 recon_y /= total_E;
181 recon_phi = atan2(recon_y,recon_x);
182 recon_theta /= total_E;
183 recon_E /= total_E;
184 recon_t /= total_E;
185 }
186 }
187 // Output new shower parameters for use in creating new showers.
188 recon_showers_phi.push_back(recon_phi);
189 recon_showers_theta.push_back(recon_theta*180.0/TMath::Pi());
190 recon_showers_E.push_back(recon_E);
191 recon_showers_t.push_back(recon_t);
192 bcal_showers.erase(bcal_showers.begin());
193 }
194
195//Create new showers. Grab curvature parameters based on theta and energy of the old showers.
196 for (int m = 0; m < (int)recon_showers_phi.size(); m++){ // Loop over our newly output shower parameters.
197 // First, figure out which theta and energy bins we need from the tables.
198 // We have to be a little careful with the first and last bins in both
199 // theta and energy, since we want to interpolate between two adjacent
200 // bin values. The following carefully sets the two adjacent bins based
201 // on the theta and energy values of the IU shower.
202 // We do an interpolation between values in adjacent theta bins, then
203 // a similar interpolation in adjacent energy bins so that our final z-position
204 // will be a more-or-less smooth function of both theta and energy.
205 temptheta = recon_showers_theta[m];
206 tempenergy = recon_showers_E[m];
207 k = l = 1;
208 bin = 0;
209 while (temptheta > 15){
210 temptheta -= 10;
211 bin++;
212 }
213 k = bin;
214 if (k > 11) k = 11;
215 bin = 0;
216 while (tempenergy > 75){
217 tempenergy -= 50;
218 bin++;
219 }
220 l = bin;
221 if (l > 31) l = 31;
222
223 if (k == 11){
224 k2 = k;
225 k--;
226 }
227 else if ((recon_showers_theta[m] - 10*(k+1)) > 0 || k == 0) k2 = k + 1;
228 else {
229 k2 = k;
230 k--;
231 }
232 if (l == 31){
233 l2 = l;
234 l--;
235 }
236 else if ((recon_showers_E[m] - 50*(l+1)) > 0 || l == 0) l2 = l + 1;
237 else {
238 l2 = l;
239 l--;
240 }
241
242 // Now that we have the bins, we can use the position and sigma arrays to extrapolate or interpolate from the nearest two bins.
243 // We have to again be a bit careful with the first and last angle bins, since they don't follow the regular 10 degree
244 // interval pattern, sitting at 5 and 115 degrees respectively. We find two interpolation offsets for the position and two
245 // for the sigma for each layer. These are then added to the positions and widths to give us our final z-bin position with error
246 // for each layer and both central and side sectors (8 pairs in all).
247 for (int ii = 0; ii < 2; ii++){
248 for (int jj = 0; jj < 4; jj++){
249 if (k == 0){ // Treat the k = 0 bin differently. It lies at 15 degrees rather than the expected (from the pattern) 10 degrees.
250 posoffset1 = (position[ii][jj][k][l] - position[ii][jj][k2][l])/(-5)*(recon_showers_theta[m] - 15);
251 sigoffset1 = (sigma[ii][jj][k][l] - sigma[ii][jj][k2][l])/(-5)*(recon_showers_theta[m] - 15);
252 }
253 else if (k2 == 11){ // Treat the k = 11 bin differently. It lies at 115 degrees rather than the expected (from the pattern) 120 degrees.
254 posoffset1 = (position[ii][jj][k][l] - position[ii][jj][k2][l])/(-5)*(recon_showers_theta[m] - 110);
255 sigoffset1 = (sigma[ii][jj][k][l] - sigma[ii][jj][k2][l])/(-5)*(recon_showers_theta[m] - 110);
256 }
257 else {
258 posoffset1 = (position[ii][jj][k][l] - position[ii][jj][k2][l])/(-10)*(recon_showers_theta[m] - 10*((double)k+1));
259 sigoffset1 = (sigma[ii][jj][k][l] - sigma[ii][jj][k2][l])/(-10)*(recon_showers_theta[m] - 10*((double)k+1));
260 }
261 posoffset2 = (position[ii][jj][k][l] - position[ii][jj][k][l2])/(-50)*(recon_showers_E[m] - 50*((double)l+1));
262 sigoffset2 = (sigma[ii][jj][k][l] - sigma[ii][jj][k][l2])/(-50)*(recon_showers_E[m] - 50*((double)l+1));
263 zbinposition[ii][jj] = position[ii][jj][k][l] + posoffset1 + posoffset2;
264 zbinwidth[ii][jj] = (sigma[ii][jj][k][l] + sigoffset1 + sigoffset2);
265 }
266 }
267
268 // We can now create the new shower.
269 DBCALShower *shower = new DBCALShower;
270
271 for (int ii = 0; ii < (int)Points.size(); ii++){
272 double delPhi;
273 if (Points[ii]->phi() > TMath::Pi()) delPhi = recon_showers_phi[m] - (Points[ii]->phi() - 2*TMath::Pi()); // Points->Phi() is measured from 0 to 2pi rather than from -pi to pi.
274 else delPhi = recon_showers_phi[m] - Points[ii]->phi();
275 if (delPhi > TMath::Pi()) delPhi -= 2*TMath::Pi();
276 if (delPhi < -1*TMath::Pi()) delPhi += 2*TMath::Pi();
277 if (fabs(delPhi) > PHITHRESHOLD) continue; // Don't include the point in this shower if it is far away in phi.
278 if (fabs(delPhi) < 0.020452475) i = 0; // 5/8 of a column width in phi. This is a 'central' shower (i = 1).
279 else i = 1;
280 j = Points[ii]->layer() - 1;
281
282 // Below, we compare the point to the shower to see if it should be included in the shower.
283 // This uses the thresholds defined at the top of this document. However, we have to be careful
284 // around the ends of the calorimeter. Points that are reconstructed outside of the calorimeter
285 // need to still be considered. So, for theta >= 119 and for theta <= 13, we must include all
286 // points past the end of the BCAL. Also, for theta <= 17, we increase the z-bin width by a factor
287 // of 1.7. These angles and widths might have to be played with a bit to find the best values,
288 // but for now, this seems to give good results (meaning, we catch most of the points in the
289 // event and put them in reasonable showers).
290 if (recon_showers_theta[m] >= 119.){ // Get all points past the edge of the calorimeter for showers near the edge.
291 if (((fabs(Points[ii]->z() - zbinposition[i][j]) < (zbinwidth[i][j] + ZTHRESHOLD)) || (Points[ii]->z() < -48)) && (fabs(Points[ii]->t() - recon_showers_t[m]) < TTHRESHOLD) && (Points[ii]->E() > ETHRESHOLD)){
292 CurvaturePoints.push_back(Points[ii]);
293 overlap.push_back(ii);
294 }
295 }
296 else if (recon_showers_theta[m] <= 13.){ // Get all points past the edge of the calorimeter for showers near the edge.
297 if (((fabs(Points[ii]->z() - zbinposition[i][j]) < 1.7*(zbinwidth[i][j] + ZTHRESHOLD)) || (Points[ii]->z() > 342)) && (fabs(Points[ii]->t() - recon_showers_t[m]) < TTHRESHOLD) && (Points[ii]->E() > ETHRESHOLD)){
298 CurvaturePoints.push_back(Points[ii]);
299 overlap.push_back(ii);
300 }
301 }
302 else if (recon_showers_theta[m] <= 17.){ // Use a larger z-bin-width for more forward-directed showers.
303 if ((fabs(Points[ii]->z() - zbinposition[i][j]) < 1.7*(zbinwidth[i][j] + ZTHRESHOLD)) && (fabs(Points[ii]->t() - recon_showers_t[m]) < TTHRESHOLD) && (Points[ii]->E() > ETHRESHOLD)){
304 CurvaturePoints.push_back(Points[ii]);
305 overlap.push_back(ii);
306 }
307 }
308 else if (recon_showers_theta[m] <= 25.){ // Use a larger z-bin-width for more forward-directed showers.
309 if ((fabs(Points[ii]->z() - zbinposition[i][j]) < 1.4*(zbinwidth[i][j] + ZTHRESHOLD)) && (fabs(Points[ii]->t() - recon_showers_t[m]) < TTHRESHOLD) && (Points[ii]->E() > ETHRESHOLD)){
310 CurvaturePoints.push_back(Points[ii]);
311 overlap.push_back(ii);
312 }
313 }
314 else if (recon_showers_theta[m] <= 50.){ // Use a larger z-bin-width for more forward-directed showers.
315 if ((fabs(Points[ii]->z() - zbinposition[i][j]) < 1.2*(zbinwidth[i][j] + ZTHRESHOLD)) && (fabs(Points[ii]->t() - recon_showers_t[m]) < TTHRESHOLD) && (Points[ii]->E() > ETHRESHOLD)){
316 CurvaturePoints.push_back(Points[ii]);
317 overlap.push_back(ii);
318 }
319 }
320 else if ((fabs(Points[ii]->z() - zbinposition[i][j]) < (zbinwidth[i][j] + ZTHRESHOLD)) && (fabs(Points[ii]->t() - recon_showers_t[m]) < TTHRESHOLD) && (Points[ii]->E() > ETHRESHOLD)){
321 CurvaturePoints.push_back(Points[ii]);
322 overlap.push_back(ii);
323 }
324 }
325
326 // Average out the showers. This is done in much the same way as the KLOE code was.
327 // x, y, z, and t are averaged weighted by the energy of the point squared, and
328 // sig_x, sig_y, sig_z, and sig_t are calculated alongside them.
329 double x=0,y=0,z=0,t=0;
330 int N_cell=0;
331 double sig_x=0,sig_y=0,sig_z=0,sig_t=0;
332 double total_E = 0; // Total Energy.
333 double total_E2 = 0; // Total Energy squared. Weight averages by E^2.
334 double total_E2_squared = 0; // Used for calculating n_eff.
335 double wt = 0; // Save writing: wt = CurvaturePoints[ii]->E().
336 bool average_layer4 = false; // Keep track of showers with only layer four points.
337 int n = CurvaturePoints.size();
338 int n4 = 0; // Number of layer four points in the shower.
339
340 for (int ii = 0; ii < (int)CurvaturePoints.size(); ii++){
341 if (CurvaturePoints[ii]->layer() == 4) n4++;
342 }
343 if (n == n4) average_layer4 = true; // If all points are layer four points, include layer four points in the averages below.
344
345 for (int ii = 0; ii < (int)CurvaturePoints.size(); ii++){
346 if (CurvaturePoints[ii]->layer() != 4 || average_layer4) wt = CurvaturePoints[ii]->E();
347 else wt = 0;
348 total_E += CurvaturePoints[ii]->E(); // We still count energy from points in layer four.
349 total_E2 += wt*wt;
350 total_E2_squared += wt*wt*wt*wt;
351 x += CurvaturePoints[ii]->r()*cos(CurvaturePoints[ii]->phi())*wt*wt;
352 y += CurvaturePoints[ii]->r()*sin(CurvaturePoints[ii]->phi())*wt*wt;
353 sig_x += CurvaturePoints[ii]->r()*CurvaturePoints[ii]->r()*cos(CurvaturePoints[ii]->phi())*cos(CurvaturePoints[ii]->phi())*wt*wt;
354 sig_y += CurvaturePoints[ii]->r()*CurvaturePoints[ii]->r()*sin(CurvaturePoints[ii]->phi())*sin(CurvaturePoints[ii]->phi())*wt*wt;
355 z += CurvaturePoints[ii]->z()*wt*wt;
356 sig_z += CurvaturePoints[ii]->z()*CurvaturePoints[ii]->z()*wt*wt;
357 t += CurvaturePoints[ii]->t()*wt*wt;
358 sig_t += CurvaturePoints[ii]->t()*CurvaturePoints[ii]->t()*wt*wt;
359 N_cell++;
360 shower->AddAssociatedObject(CurvaturePoints[ii]);
361 }
362 while (overlap.size() != 0){
363 Points.erase(Points.begin()+overlap.back());
364 overlap.pop_back();
365 }
366 CurvaturePoints.clear();
367
368 // the variance of the mean of a weighted distribution is s^2/n_eff, where s^2 is the variance of the sample and n_eff is as calculated below.
369 double n_eff = total_E2*total_E2/total_E2_squared;
370
371 x /= total_E2;
372 sig_x /= total_E2;
373 sig_x = sqrt(sig_x - x*x)/sqrt(n_eff);
374 y /= total_E2;
375 sig_y /= total_E2;
376 sig_y = sqrt(sig_y - y*y)/sqrt(n_eff);
Value stored to 'sig_y' is never read
377 z /= total_E2;
378 sig_z /= total_E2;
379 sig_z = sqrt(sig_z - z*z)/sqrt(n_eff);
380 t /= total_E2;
381 sig_t /= total_E2;
382 sig_t = sqrt(sig_t - t*t)/sqrt(n_eff);
383
384 // Force showers to be inside the BCal's z-coordinate range.
385 double bcal_down = DBCALGeometry::GetBCAL_center() + DBCALGeometry::GetBCAL_length()/2.0 - m_zTarget;
386 double bcal_up = DBCALGeometry::GetBCAL_center() - DBCALGeometry::GetBCAL_length()/2.0 - m_zTarget;
387 if (z > bcal_down) z = bcal_down;
388 if (z < bcal_up) z = bcal_up;
389
390 shower->E_raw = total_E;
391 shower->x = x;
392 shower->y = y;
393 shower->z = z + m_zTarget;
394 shower->t = t;
395 shower->N_cell = N_cell;
396 // shower->xErr = sig_x;
397 // shower->yErr = sig_y;
398 // shower->zErr = sig_z;
399 // shower->tErr = sig_t;
400
401 // calibrate energy:
402 // Energy calibration has a z dependence -- the
403 // calibration comes from fitting E_rec / E_gen to scale * E_gen^nonlin
404 // for slices of z. These fit parameters (scale and nonlin) are then plotted
405 // as a function of z and fit.
406 float r = sqrt( shower->x * shower->x + shower->y * shower->y );
407 float zEntry = ( shower->z - m_zTarget ) * ( DBCALGeometry::GetBCAL_inner_rad() / r );
408 float scale = m_scaleZ_p0 + m_scaleZ_p1*zEntry + m_scaleZ_p2*(zEntry*zEntry) + m_scaleZ_p3*(zEntry*zEntry*zEntry);
409 float nonlin = m_nonlinZ_p0 + m_nonlinZ_p1*zEntry + m_nonlinZ_p2*(zEntry*zEntry) + m_nonlinZ_p3*(zEntry*zEntry*zEntry);
410
411 shower->E = pow( (shower->E_raw ) / scale, 1 / ( 1 + nonlin ) );
412
413 //copy xyz errors into covariance matrix
414 // shower->xyzCovariance.ResizeTo(3,3);
415 // shower->xyzCovariance[0][0] = shower->xErr*shower->xErr;
416 // shower->xyzCovariance[1][1] = shower->yErr*shower->yErr;
417 // shower->xyzCovariance[2][2] = shower->zErr*shower->zErr;
418
419 _data.push_back(shower);
420 }
421
422 return NOERROR;
423}
424
425//------------------
426// erun
427//------------------
428jerror_t DBCALShower_factory_CURVATURE::erun(void)
429{
430 return NOERROR;
431}
432
433//------------------
434// fini
435//------------------
436jerror_t DBCALShower_factory_CURVATURE::fini(void)
437{
438 return NOERROR;
439}
440