1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | |
14 | |
15 | |
16 | |
17 | |
18 | |
19 | |
20 | |
21 | |
22 | |
23 | |
24 | |
25 | |
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 | |
43 | |
44 | jerror_t DBCALShower_factory_CURVATURE::init(void) |
45 | { |
46 | if( ! DBCALGeometry::summingOn() ) { |
47 | |
48 | |
49 | |
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 | |
64 | |
65 | |
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 | |
82 | |
83 | jerror_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 | |
92 | |
93 | |
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){ |
108 | position[ii][layer - 1][11][energy/50 - 1] = dataposition - 48; |
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; |
113 | sigma[ii][layer - 1][angle/10 - 1][energy/50 - 1] = datasigma; |
114 | } |
115 | } |
116 | curvature_parameters.clear(); |
117 | } |
118 | |
119 | |
120 | PHITHRESHOLD = 4*0.06544792; |
121 | ZTHRESHOLD = 35.*k_cm; |
122 | TTHRESHOLD = 8.0*k_nsec; |
123 | ETHRESHOLD = 1.0*k_keV; |
124 | |
125 | return NOERROR; |
126 | } |
127 | |
128 | |
129 | |
130 | |
131 | jerror_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 |
|
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; |
152 | double total_E = bcal_showers[0]->E; |
153 | for (int ii = 1; ii < (int)bcal_showers.size(); ii++){ |
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); |
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){ |
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; |
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)*bcal_showers[overlap.back()]->E; |
174 | total_E += bcal_showers[overlap.back()]->E; |
175 | |
176 | bcal_showers.erase(bcal_showers.begin()+overlap.back()); |
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 | |
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 | |
196 | for (int m = 0; m < (int)recon_showers_phi.size(); m++){ |
197 | |
198 | |
199 | |
200 | |
201 | |
202 | |
203 | |
204 | |
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 | |
243 | |
244 | |
245 | |
246 | |
247 | for (int ii = 0; ii < 2; ii++){ |
248 | for (int jj = 0; jj < 4; jj++){ |
249 | if (k == 0){ |
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){ |
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 | |
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()); |
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; |
278 | if (fabs(delPhi) < 0.020452475) i = 0; |
279 | else i = 1; |
280 | j = Points[ii]->layer() - 1; |
281 | |
282 | |
283 | |
284 | |
285 | |
286 | |
287 | |
288 | |
289 | |
290 | if (recon_showers_theta[m] >= 119.){ |
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.){ |
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.){ |
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.){ |
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.){ |
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 | |
327 | |
328 | |
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; |
333 | double total_E2 = 0; |
334 | double total_E2_squared = 0; |
335 | double wt = 0; |
336 | bool average_layer4 = false; |
337 | int n = CurvaturePoints.size(); |
338 | int n4 = 0; |
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; |
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(); |
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 | |
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); |
377 | z /= total_E2; |
378 | sig_z /= total_E2; |
379 | sig_z = sqrt(sig_z - z*z)/sqrt(n_eff); |
| Value stored to 'sig_z' is never read |
380 | t /= total_E2; |
381 | sig_t /= total_E2; |
382 | sig_t = sqrt(sig_t - t*t)/sqrt(n_eff); |
383 | |
384 | |
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 | |
397 | |
398 | |
399 | |
400 | |
401 | |
402 | |
403 | |
404 | |
405 | |
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 | |
414 | |
415 | |
416 | |
417 | |
418 | |
419 | _data.push_back(shower); |
420 | } |
421 | |
422 | return NOERROR; |
423 | } |
424 | |
425 | |
426 | |
427 | |
428 | jerror_t DBCALShower_factory_CURVATURE::erun(void) |
429 | { |
430 | return NOERROR; |
431 | } |
432 | |
433 | |
434 | |
435 | |
436 | jerror_t DBCALShower_factory_CURVATURE::fini(void) |
437 | { |
438 | return NOERROR; |
439 | } |
440 | |