3 #include <phparameter/PHParameterInterface.h>
40 for (
int j = 0; j <
nRadii; j++)
57 CalculateVertices(
nStripes_R1,
nPads_R1,
R1_e,
spacing_R1_e,
x1a_R1_e,
y1a_R1_e,
x1b_R1_e,
y1b_R1_e,
x2a_R1_e,
y2a_R1_e,
x2b_R1_e,
y2b_R1_e,
x3a_R1_e,
y3a_R1_e,
x3b_R1_e,
y3b_R1_e,
padfrac_R1,
str_width_R1_e,
widthmod_R1_e,
nGoodStripes_R1_e,
keepUntil_R1_e,
nStripesIn_R1_e,
nStripesBefore_R1_e);
58 CalculateVertices(
nStripes_R1,
nPads_R1,
R1,
spacing_R1,
x1a_R1,
y1a_R1,
x1b_R1,
y1b_R1,
x2a_R1,
y2a_R1,
x2b_R1,
y2b_R1,
x3a_R1,
y3a_R1,
x3b_R1,
y3b_R1,
padfrac_R1,
str_width_R1,
widthmod_R1,
nGoodStripes_R1,
keepUntil_R1,
nStripesIn_R1,
nStripesBefore_R1);
59 CalculateVertices(
nStripes_R2,
nPads_R2,
R2,
spacing_R2,
x1a_R2,
y1a_R2,
x1b_R2,
y1b_R2,
x2a_R2,
y2a_R2,
x2b_R2,
y2b_R2,
x3a_R2,
y3a_R2,
x3b_R2,
y3b_R2,
padfrac_R2,
str_width_R2,
widthmod_R2,
nGoodStripes_R2,
keepUntil_R2,
nStripesIn_R2,
nStripesBefore_R2);
60 CalculateVertices(
nStripes_R3,
nPads_R3,
R3,
spacing_R3,
x1a_R3,
y1a_R3,
x1b_R3,
y1b_R3,
x2a_R3,
y2a_R3,
x2b_R3,
y2b_R3,
x3a_R3,
y3a_R3,
x3b_R3,
y3b_R3,
padfrac_R3,
str_width_R3,
widthmod_R3,
nGoodStripes_R3,
keepUntil_R3,
nStripesIn_R3,
nStripesBefore_R3);
88 std::cout <<
"PHG4TpcCentralMembrane::InitRun - electrons_per_stripe: " <<
electrons_per_stripe << std::endl;
89 std::cout <<
"PHG4TpcCentralMembrane::InitRun - electrons_per_gev " <<
electrons_per_gev << std::endl;
93 auto* g4hit = findNode::getClass<PHG4HitContainer>(topNode,
hitnodename.c_str());
96 std::cout <<
Name() <<
" Could not locate G4HIT node " <<
hitnodename << std::endl;
119 for (
int i = 0; i < 18; i++)
121 for (
int j = 0; j < 8; j++)
151 for (
const auto& hit : PHG4Hits)
168 auto g4hitcontainer = findNode::getClass<PHG4HitContainer>(topNode,
hitnodename.c_str());
179 g4hitcontainer->AddHit(detId,
copy);
194 static constexpr
double Ne_dEdx = 1.56;
195 static constexpr
double CF4_dEdx = 7.00;
196 static constexpr
double Ne_NTotal = 43;
197 static constexpr
double CF4_NTotal = 100;
198 static constexpr
double Tpc_NTot = 0.5 * Ne_NTotal + 0.5 * CF4_NTotal;
199 static constexpr
double Tpc_dEdx = 0.5 * Ne_dEdx + 0.5 * CF4_dEdx;
200 static constexpr
double Tpc_ElectronsPerKeV = Tpc_NTot / Tpc_dEdx;
211 int nStripes,
int nPads,
212 const std::array<double, nRadii>&
R,
213 std::array<double, nRadii>& spacing,
214 double x1a[][nRadii],
double y1a[][nRadii],
215 double x1b[][nRadii],
double y1b[][nRadii],
216 double x2a[][nRadii],
double y2a[][nRadii],
217 double x2b[][nRadii],
double y2b[][nRadii],
218 double x3a[][nRadii],
double y3a[][nRadii],
219 double x3b[][nRadii],
double y3b[][nRadii],
221 double str_width[][nRadii],
222 const std::array<double, nRadii>& widthmod,
223 std::array<int, nRadii>& nGoodStripes,
224 const std::array<int, nRadii>& keepUntil,
225 std::array<int, nRadii>& nStripesIn,
226 std::array<int, nRadii>& nStripesBefore)
228 const double phi_module =
M_PI / 6.0;
229 const int pr_mult = 3;
230 const int dw_mult = 8;
231 const double diffwidth = 0.6 *
mm;
232 const double adjust = 0.015;
248 for (
int i = 0; i <
nRadii; i++)
250 spacing[i] = 2.0 * ((dw_mult * diffwidth / R[i]) + (pr_mult * phi_module / nPads));
254 for (
int j = 0; j <
nRadii; j++)
261 theta = i * spacing[j] + (spacing[j] / 2) - adjust;
262 cx[i_out][j] = R[j] * cos(theta);
263 cy[i_out][j] = R[j] * sin(theta);
267 theta = (i + 1) * spacing[j] - adjust;
268 cx[i_out][j] = R[j] * cos(theta);
269 cy[i_out][j] = R[j] * sin(theta);
273 corner[0].SetXYZ(-padfrac +
arc_r, -(widthmod[j] * str_width[i][j]) / 2, 0);
274 corner[1].SetXYZ(padfrac -
arc_r, -(widthmod[j] * str_width[i][j]) / 2, 0);
275 corner[2].SetXYZ(-padfrac +
arc_r, (widthmod[j] * str_width[i][j]) / 2, 0);
276 corner[3].SetXYZ(padfrac -
arc_r, (widthmod[j] * str_width[i][j]) / 2, 0);
278 TVector3 rotatedcorner[4];
279 for (
int i = 0; i < 4; i++)
281 rotatedcorner[i] = corner[i];
282 rotatedcorner[i].RotateZ(theta);
285 x1a[i_out][j] = rotatedcorner[0].X() + cx[i_out][j];
286 x1b[i_out][j] = rotatedcorner[1].X() + cx[i_out][j];
287 x2a[i_out][j] = rotatedcorner[2].X() + cx[i_out][j];
288 x2b[i_out][j] = rotatedcorner[3].X() + cx[i_out][j];
290 y1a[i_out][j] = rotatedcorner[0].Y() + cy[i_out][j];
291 y1b[i_out][j] = rotatedcorner[1].Y() + cy[i_out][j];
292 y2a[i_out][j] = rotatedcorner[2].Y() + cy[i_out][j];
293 y2b[i_out][j] = rotatedcorner[3].Y() + cy[i_out][j];
331 x3a[i_out][j] = (x1a[i_out][j] + x2a[i_out][j]) / 2.0;
332 y3a[i_out][j] = (y1a[i_out][j] + y2a[i_out][j]) / 2.0;
333 x3b[i_out][j] = (x1b[i_out][j] + x2b[i_out][j]) / 2.0;
334 y3b[i_out][j] = (y1b[i_out][j] + y2b[i_out][j]) / 2.0;
342 nStripesBefore[j] = nStripesIn[j - 1] + nStripesBefore[j - 1];
345 nGoodStripes[j] = i_out;
352 TVector3 dummyPos0, dummyPos1;
365 else if (moduleID == 1)
372 else if (moduleID == 2)
379 else if (moduleID == 3)
386 botvert->
set_z(0, 0.0);
394 TVector3 dummyPos0, dummyPos1;
407 else if (moduleID == 1)
414 else if (moduleID == 2)
421 else if (moduleID == 3)
428 topvert->
set_z(0, 0.0);
434 const double x1a[][nRadii],
const double x1b[][nRadii],
435 const double x2a[][nRadii],
const double x2b[][nRadii],
436 const double y1a[][nRadii],
const double y1b[][nRadii],
437 const double y2a[][nRadii],
const double y2b[][nRadii],
438 const double x3a[][nRadii],
const double y3a[][nRadii],
439 const double x3b[][nRadii],
const double y3b[][nRadii],
440 double x,
double y,
const std::array<int, nRadii>& nGoodStripes)
const
444 for (
int j = 0; j <
nRadii; j++)
446 for (
int i = 0; i < nGoodStripes[j]; i++)
448 if (((y1a[i][j] > y) != (y2a[i][j] > y) && (x < (x2a[i][j] - x1a[i][j]) * (y - y1a[i][j]) / (y2a[i][j] - y1a[i][j]) + x1a[i][j])))
450 if (((y1b[i][j] > y) != (y1a[i][j] > y) && (x < (x1a[i][j] - x1b[i][j]) * (y - y1b[i][j]) / (y1a[i][j] - y1b[i][j]) + x1b[i][j])))
452 if (((y2b[i][j] > y) != (y1b[i][j] > y) && (x < (x1b[i][j] - x2b[i][j]) * (y - y2b[i][j]) / (y1b[i][j] - y2b[i][j]) + x2b[i][j])))
454 if (((y2a[i][j] > y) != (y2b[i][j] > y) && (x < (x2b[i][j] - x2a[i][j]) * (y - y2a[i][j]) / (y2b[i][j] - y2a[i][j]) + x2a[i][j])))
460 if (((x - x3a[i][j]) * (x - x3a[i][j]) + (y - y3a[i][j]) * (y - y3a[i][j])) <=
arc_r *
arc_r)
464 else if (((x - x3b[i][j]) * (x - x3b[i][j]) + (y - y3b[i][j]) * (y - y3b[i][j])) <= arc_r * arc_r)
476 const double phi_petal =
M_PI / 9.0;
477 const double end_R1_e = 312.0 *
mm;
478 const double end_R1 = 408.0 *
mm;
479 const double end_R2 = 580.0 *
mm;
481 double r,
phi, phimod, xmod, ymod;
483 r = sqrt(xcheck * xcheck + ycheck * ycheck);
484 phi = atan(ycheck / xcheck);
485 if ((xcheck < 0.0) && (ycheck > 0.0))
489 else if ((xcheck > 0.0) && (ycheck < 0.0))
491 phi = phi + 2.0 *
M_PI;
494 phimod = fmod(phi, phi_petal);
495 xmod = r * cos(phimod);
496 ymod = r * sin(phimod);
502 result =
SearchModule(
nStripes_R1,
x1a_R1_e,
x1b_R1_e,
x2a_R1_e,
x2b_R1_e,
y1a_R1_e,
y1b_R1_e,
y2a_R1_e,
y2b_R1_e,
x3a_R1_e,
y3a_R1_e,
x3b_R1_e,
y3b_R1_e, xmod, ymod,
nGoodStripes_R1_e);
504 else if ((r > end_R1_e) && (r <= end_R1))
506 result =
SearchModule(
nStripes_R1,
x1a_R1,
x1b_R1,
x2a_R1,
x2b_R1,
y1a_R1,
y1b_R1,
y2a_R1,
y2b_R1,
x3a_R1,
y3a_R1,
x3b_R1,
y3b_R1, xmod, ymod,
nGoodStripes_R1);
508 else if ((r > end_R1) && (r <= end_R2))
510 result =
SearchModule(
nStripes_R2,
x1a_R2,
x1b_R2,
x2a_R2,
x2b_R2,
y1a_R2,
y1b_R2,
y2a_R2,
y2b_R2,
x3a_R2,
y3a_R2,
x3b_R2,
y3b_R2, xmod, ymod,
nGoodStripes_R2);
512 else if ((r > end_R2) && (r <=
end_CM))
514 result =
SearchModule(
nStripes_R3,
x1a_R3,
x1b_R3,
x2a_R3,
x2b_R3,
y1a_R3,
y1b_R3,
y2a_R3,
y2b_R3,
x3a_R3,
y3a_R3,
x3b_R3,
y3b_R3, xmod, ymod,
nGoodStripes_R3);
522 const double phi_petal =
M_PI / 9.0;
524 TVector3 dummyPos0, dummyPos1;
539 else if (moduleID == 1)
544 else if (moduleID == 2)
549 else if (moduleID == 3)
560 dummyPos0.RotateZ(petalID * phi_petal);
561 hit->
set_x(0, dummyPos0.X());
562 hit->
set_y(0, dummyPos0.Y());
584 else if (moduleID == 1)
589 else if (moduleID == 2)
594 else if (moduleID == 3)
605 dummyPos1.RotateZ(petalID * phi_petal);
606 hit->
set_x(1, dummyPos1.X());
607 hit->
set_y(1, dummyPos1.Y());
653 int result, rID, petalID;
659 const double phi_petal =
M_PI / 9.0;
661 double r,
phi, phimod, xmod, ymod;
671 r = sqrt(xcheck * xcheck + ycheck * ycheck);
672 phi = atan(ycheck / xcheck);
673 if ((xcheck < 0.0) && (ycheck > 0.0))
677 else if ((xcheck > 0.0) && (ycheck < 0.0))
679 phi = phi + 2.0 *
M_PI;
682 phimod = fmod(phi, phi_petal);
683 xmod = r * cos(phimod);
684 ymod = r * sin(phimod);
686 petalID = phi / phi_petal;
688 for (
int j = 0; j <
nRadii; j++)
693 std::cout <<
"rID: " << rID << std::endl;
713 dist = fabs((-m) * xmod + ymod) / sqrt(1 + m * m);
734 dist = fabs(m * xmod - ymod) / sqrt(1 + m * m);
751 dist = fabs(m * xmod - ymod) / sqrt(1 + m * m);
768 dist = fabs(m * xmod - ymod) / sqrt(1 + m * m);