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;
54 m_nonlinZ_p0 = -0.0147086;
55 m_nonlinZ_p1 = 9.69207e-05;
71 vector< vector<double> > curvature_parameters;
76 for (
int ii = 0; ii < 2; ii++){
78 if (loop->GetCalib(
"/BCAL/curvature_central", curvature_parameters)) jout <<
"Error loading /BCAL/curvature_central !" << endl;
81 if (loop->GetCalib(
"/BCAL/curvature_side", curvature_parameters)) jout <<
"Error loading /BCAL/curvature_side !" << endl;
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];
90 position[ii][
layer - 1][11][energy/50 - 1] = dataposition - 48;
91 sigma[ii][
layer - 1][11][energy/50 - 1] = datasigma;
94 position[ii][
layer - 1][angle/10 - 1][energy/50 - 1] = dataposition - 48;
95 sigma[ii][
layer - 1][angle/10 - 1][energy/50 - 1] = datasigma;
98 curvature_parameters.clear();
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];
109 PHITHRESHOLD = 4*0.06544792;
110 ZTHRESHOLD = 35.*
k_cm;
112 ETHRESHOLD = 1.0*
k_keV;
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;
130 loop->Get(bcal_showers,
"IU");
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;
141 double total_E = bcal_showers[0]->E;
142 for (
int ii = 1; ii < (int)bcal_showers.size(); ii++){
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);
150 if (overlap.size() != 0){
153 recon_theta *= recon_E;
157 while (overlap.size() != 0){
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;
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)*bcal_showers[overlap.back()]->E;
163 total_E += bcal_showers[overlap.back()]->E;
165 bcal_showers.erase(bcal_showers.begin()+overlap.back());
167 if (overlap.size() == 0){
170 recon_phi = atan2(recon_y,recon_x);
171 recon_theta /= total_E;
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());
185 for (
int m = 0; m < (int)recon_showers_phi.size(); m++){
194 temptheta = recon_showers_theta[m];
195 tempenergy = recon_showers_E[m];
198 while (temptheta > 15){
205 while (tempenergy > 75){
216 else if ((recon_showers_theta[m] - 10*(k+1)) > 0 || k == 0) k2 = k + 1;
225 else if ((recon_showers_E[m] - 50*(l+1)) > 0 || l == 0) l2 = l + 1;
236 for (
int ii = 0; ii < 2; ii++){
237 for (
int jj = 0; jj < 4; jj++){
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);
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);
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));
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);
260 for (
int ii = 0; ii < (int)Points.size(); ii++){
262 if (Points[ii]->phi() > TMath::Pi()) delPhi = recon_showers_phi[m] - (Points[ii]->phi() - 2*TMath::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;
267 if (fabs(delPhi) < 0.020452475) i = 0;
269 j = Points[ii]->layer() - 1;
279 if (recon_showers_theta[m] >= 119.){
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);
285 else if (recon_showers_theta[m] <= 13.){
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);
291 else if (recon_showers_theta[m] <= 17.){
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);
297 else if (recon_showers_theta[m] <= 25.){
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);
303 else if (recon_showers_theta[m] <= 50.){
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);
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);
318 double x=0,
y=0,z=0,t=0;
320 double sig_x=0,sig_y=0,sig_z=0,sig_t=0;
323 double total_E2_squared = 0;
325 bool average_layer4 =
false;
326 int n = CurvaturePoints.size();
329 for (
int ii = 0; ii < (int)CurvaturePoints.size(); ii++){
330 if (CurvaturePoints[ii]->
layer() == 4) n4++;
332 if (n == n4) average_layer4 =
true;
334 for (
int ii = 0; ii < (int)CurvaturePoints.size(); ii++){
335 if (CurvaturePoints[ii]->
layer() != 4 || average_layer4) wt = CurvaturePoints[ii]->E();
337 total_E += CurvaturePoints[ii]->E();
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;
349 shower->AddAssociatedObject(CurvaturePoints[ii]);
351 while (overlap.size() != 0){
352 Points.erase(Points.begin()+overlap.back());
355 CurvaturePoints.clear();
358 double n_eff = total_E2*total_E2/total_E2_squared;
362 sig_x =
sqrt(sig_x - x*x)/
sqrt(n_eff);
368 sig_z =
sqrt(sig_z - z*z)/
sqrt(n_eff);
371 sig_t =
sqrt(sig_t - t*t)/
sqrt(n_eff);
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;
379 shower->
E_raw = total_E;
382 shower->
z = z + m_zTarget;
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);
400 shower->
E = pow( (shower->
E_raw ) / scale, 1 / ( 1 + nonlin ) );
408 _data.push_back(shower);
jerror_t init(void)
Called once at program start.
DGeometry * GetDGeometry(unsigned int run_number)
Double_t sigma[NCHANNELS]
jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
jerror_t erun(void)
Called everytime run number changes, provided brun has been called.
bool GetTargetZ(double &z_target) const
z-location of center of target
jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
jerror_t fini(void)
Called after last event of last event source has been processed.