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