File indexing completed on 2025-08-05 08:17:52
0001 #include "PHG4TpcCylinderGeom.h"
0002 #include "PHG4CylinderCellDefs.h"
0003
0004 #include <phool/phool.h>
0005
0006 #include <cstdlib>
0007
0008 namespace
0009 {
0010
0011 using array_t = std::array<std::vector<double>, PHG4TpcCylinderGeom::NSides>;
0012 std::ostream& operator<<(std::ostream& out, const array_t& array)
0013 {
0014 out << "{ ";
0015 for (const auto& iside : array)
0016 {
0017 out << "{";
0018 bool first = true;
0019 for (const auto& value : iside)
0020 {
0021 if (!first)
0022 {
0023 out << ", ";
0024 }
0025 first = false;
0026 out << value;
0027 }
0028 out << "} ";
0029 }
0030 out << " }";
0031 return out;
0032 }
0033 }
0034
0035 std::ostream& operator<<(std::ostream& out, const PHG4TpcCylinderGeom& geom)
0036 {
0037 out << "PHG4TpcCylinderGeom - layer: " << geom.layer << std::endl;
0038 out
0039 << " binnig: " << geom.binning
0040 << ", radius: " << geom.radius
0041 << ", nzbins: " << geom.nzbins
0042 << ", zmin: " << geom.zmin
0043 << ", zstep: " << geom.zstep
0044 << ", nphibins: " << geom.nphibins
0045 << ", phimin: " << geom.phimin
0046 << ", phistep: " << geom.phistep
0047 << ", thickness: " << geom.thickness
0048 << std::endl;
0049
0050 out << " sector_R_bias: " << geom.sector_R_bias << std::endl;
0051 out << " sector_Phi_bias: " << geom.sector_Phi_bias << std::endl;
0052 out << " sector_min_Phi: " << geom.sector_min_Phi << std::endl;
0053 out << " sector_max_Phi: " << geom.sector_max_Phi << std::endl;
0054
0055 return out;
0056 }
0057
0058 void PHG4TpcCylinderGeom::set_zbins(const int i)
0059 {
0060 check_binning_method(PHG4CylinderCellDefs::sizebinning);
0061 nzbins = i;
0062 }
0063
0064 void PHG4TpcCylinderGeom::set_zmin(const double z)
0065 {
0066 check_binning_method(PHG4CylinderCellDefs::sizebinning);
0067 zmin = z;
0068 }
0069
0070 int PHG4TpcCylinderGeom::get_zbins() const
0071 {
0072 check_binning_method(PHG4CylinderCellDefs::sizebinning);
0073 return nzbins;
0074 }
0075
0076 double
0077 PHG4TpcCylinderGeom::get_zmin() const
0078 {
0079 check_binning_method(PHG4CylinderCellDefs::sizebinning);
0080 return zmin;
0081 }
0082
0083 double
0084 PHG4TpcCylinderGeom::get_zstep() const
0085 {
0086 check_binning_method(PHG4CylinderCellDefs::sizebinning);
0087 return zstep;
0088 }
0089
0090 void PHG4TpcCylinderGeom::set_zstep(const double z)
0091 {
0092 check_binning_method(PHG4CylinderCellDefs::sizebinning);
0093 zstep = z;
0094 }
0095
0096 int PHG4TpcCylinderGeom::get_phibins() const
0097 {
0098 check_binning_method_phi("PHG4TpcCylinderGeom::get_phibins");
0099 return nphibins;
0100 }
0101
0102 double
0103 PHG4TpcCylinderGeom::get_phistep() const
0104 {
0105 check_binning_method_phi("PHG4TpcCylinderGeom::get_phistep");
0106 return phistep;
0107 }
0108
0109 double
0110 PHG4TpcCylinderGeom::get_phimin() const
0111 {
0112 check_binning_method_phi("PHG4TpcCylinderGeom::get_phimin");
0113 return phimin;
0114 }
0115
0116 void PHG4TpcCylinderGeom::set_phibins(const int i)
0117 {
0118 check_binning_method_phi("PHG4TpcCylinderGeom::set_phibins");
0119 nphibins = i;
0120 }
0121
0122 void PHG4TpcCylinderGeom::set_phistep(const double r)
0123 {
0124 check_binning_method_phi("PHG4TpcCylinderGeom::set_phistep");
0125 phistep = r;
0126 }
0127
0128 void PHG4TpcCylinderGeom::set_phimin(const double r)
0129 {
0130 check_binning_method_phi("PHG4TpcCylinderGeom::set_phimin");
0131 phimin = r;
0132 }
0133
0134 int PHG4TpcCylinderGeom::get_etabins() const
0135 {
0136 check_binning_method_eta("PHG4TpcCylinderGeom::get_etabins");
0137 return nzbins;
0138 }
0139
0140 double
0141 PHG4TpcCylinderGeom::get_etastep() const
0142 {
0143 check_binning_method_eta("PHG4TpcCylinderGeom::get_etastep");
0144 return zstep;
0145 }
0146 double
0147 PHG4TpcCylinderGeom::get_etamin() const
0148 {
0149 check_binning_method_eta("PHG4TpcCylinderGeom::get_etamin");
0150 return zmin;
0151 }
0152
0153 void PHG4TpcCylinderGeom::set_etamin(const double z)
0154 {
0155 check_binning_method_eta("PHG4TpcCylinderGeom::set_etamin");
0156 zmin = z;
0157 }
0158
0159 void PHG4TpcCylinderGeom::set_etastep(const double z)
0160 {
0161 check_binning_method_eta("PHG4TpcCylinderGeom::set_etastep");
0162 zstep = z;
0163 }
0164
0165 void PHG4TpcCylinderGeom::set_etabins(const int i)
0166 {
0167 check_binning_method_eta("PHG4TpcCylinderGeom::set_etabins");
0168 nzbins = i;
0169 }
0170
0171 void PHG4TpcCylinderGeom::identify(std::ostream& os) const
0172 {
0173 os << "PHG4TpcCylinderGeom::identify - layer: " << layer
0174 << ", radius: " << radius
0175 << ", thickness: " << thickness;
0176 switch (binning)
0177 {
0178 case PHG4CylinderCellDefs::sizebinning:
0179 os << ", zbins: " << nzbins
0180 << ", zmin: " << zmin
0181 << ", zstepsize: " << zstep;
0182 break;
0183 case PHG4CylinderCellDefs::etaphibinning:
0184 case PHG4CylinderCellDefs::etaslatbinning:
0185 os << ", etabins: " << nzbins
0186 << ", etamin: " << zmin
0187 << ", etastepsize: " << zstep;
0188 break;
0189 case PHG4CylinderCellDefs::spacalbinning:
0190 os << ", etabins: " << nzbins << " for Spacal";
0191 break;
0192 default:
0193 os << "no valid binning method: " << binning << std::endl;
0194 return;
0195 break;
0196 }
0197 os << ", phimin: " << phimin
0198 << ", phibins: " << nphibins
0199 << ", phistep: " << phistep
0200 << std::endl;
0201 return;
0202 }
0203
0204 std::pair<double, double>
0205 PHG4TpcCylinderGeom::get_zbounds(const int ibin) const
0206 {
0207 if (ibin < 0 || ibin > nzbins)
0208 {
0209 std::cout << PHWHERE << " Asking for invalid bin in z: " << ibin << std::endl;
0210 exit(1);
0211 }
0212 check_binning_method(PHG4CylinderCellDefs::sizebinning);
0213 double zlow = zmin + ibin * zstep;
0214 double zhigh = zlow + zstep;
0215 return std::make_pair(zlow, zhigh);
0216 }
0217
0218 std::pair<double, double>
0219 PHG4TpcCylinderGeom::get_etabounds(const int ibin) const
0220 {
0221 if (ibin < 0 || ibin > nzbins)
0222 {
0223 std::cout << PHWHERE << " Asking for invalid bin in z: " << ibin << std::endl;
0224 exit(1);
0225 }
0226 check_binning_method_eta("PHG4TpcCylinderGeom::get_etabounds");
0227
0228 double zlow = zmin + ibin * zstep;
0229 double zhigh = zlow + zstep;
0230 return std::make_pair(zlow, zhigh);
0231 }
0232
0233 std::pair<double, double>
0234 PHG4TpcCylinderGeom::get_phibounds(const int ibin) const
0235 {
0236 if (ibin < 0 || ibin > nphibins)
0237 {
0238 std::cout << PHWHERE << "Asking for invalid bin in phi: " << ibin << std::endl;
0239 exit(1);
0240 }
0241
0242 double philow = phimin + ibin * phistep;
0243 double phihigh = philow + phistep;
0244 return std::make_pair(philow, phihigh);
0245 }
0246
0247 int PHG4TpcCylinderGeom::get_zbin(const double z) const
0248 {
0249 if (z < zmin || z > (zmin + nzbins * zstep))
0250 {
0251
0252 return -1;
0253 }
0254
0255 check_binning_method(PHG4CylinderCellDefs::sizebinning);
0256 return floor((z - zmin) / zstep);
0257 }
0258
0259 int PHG4TpcCylinderGeom::get_etabin(const double eta) const
0260 {
0261 if (eta < zmin || eta > (zmin + nzbins * zstep))
0262 {
0263
0264 return -1;
0265 }
0266 check_binning_method_eta();
0267 return floor((eta - zmin) / zstep);
0268 }
0269
0270 int PHG4TpcCylinderGeom::get_phibin_new(const double phi) const
0271 {
0272 double norm_phi = phi;
0273 if (phi < phimin || phi > (phimin + nphibins * phistep))
0274 {
0275 int nwraparound = -floor((phi - phimin) * 0.5 / M_PI);
0276 norm_phi += 2 * M_PI * nwraparound;
0277 }
0278 check_binning_method_phi();
0279 return floor((norm_phi - phimin) / phistep);
0280 }
0281
0282 int PHG4TpcCylinderGeom::find_phibin(const double phi, int side) const
0283 {
0284 double norm_phi = phi;
0285 if (phi < phimin || phi > (phimin + nphibins * phistep))
0286 {
0287 int nwraparound = -floor((phi - phimin) * 0.5 / M_PI);
0288 norm_phi += 2 * M_PI * nwraparound;
0289 }
0290
0291
0292
0293
0294
0295
0296
0297
0298 int phi_bin = -1;
0299
0300 for (std::size_t s = 0; s < sector_max_Phi[side].size(); s++)
0301 {
0302 if (norm_phi < sector_max_Phi[side][s] && norm_phi > sector_min_Phi[side][s])
0303 {
0304
0305 phi_bin = (floor(std::abs(sector_max_Phi[side][s] - norm_phi) / phistep) + nphibins / 12 * s);
0306 break;
0307 }
0308 if (s == 11)
0309 {
0310 if (norm_phi < sector_max_Phi[side][s] && norm_phi >= -M_PI)
0311 {
0312
0313 phi_bin = floor(std::abs(sector_max_Phi[side][s] - norm_phi) / phistep) + nphibins / 12 * s;
0314 break;
0315 }
0316 if (norm_phi > sector_min_Phi[side][s] + 2 * M_PI)
0317 {
0318
0319 phi_bin = floor(std::abs(sector_max_Phi[side][s] - (norm_phi - 2 * M_PI)) / phistep) + nphibins / 12 * s;
0320 break;
0321 }
0322 }
0323 }
0324 return phi_bin;
0325 }
0326
0327 float PHG4TpcCylinderGeom::get_pad_float(const double phi, int side) const
0328 {
0329 double norm_phi = phi;
0330 if (phi < phimin || phi > (phimin + nphibins * phistep))
0331 {
0332 int nwraparound = -floor((phi - phimin) * 0.5 / M_PI);
0333 norm_phi += 2 * M_PI * nwraparound;
0334 }
0335
0336
0337
0338
0339
0340
0341
0342
0343 float phi_bin = -1;
0344
0345 for (std::size_t s = 0; s < sector_max_Phi[side].size(); s++)
0346 {
0347 if (norm_phi < sector_max_Phi[side][s] && norm_phi > sector_min_Phi[side][s])
0348 {
0349
0350 phi_bin = (std::abs(sector_max_Phi[side][s] - norm_phi) / phistep) + nphibins / 12 * s;
0351 break;
0352 }
0353 if (s == 11)
0354 {
0355 if (norm_phi < sector_max_Phi[side][s] && norm_phi >= -M_PI)
0356 {
0357
0358 phi_bin = (std::abs(sector_max_Phi[side][s] - norm_phi) / phistep) + nphibins / 12 * s;
0359 break;
0360 }
0361 if (norm_phi > sector_min_Phi[side][s] + 2 * M_PI)
0362 {
0363
0364 phi_bin = (std::abs(sector_max_Phi[side][s] - (norm_phi - 2 * M_PI)) / phistep) + nphibins / 12 * s;
0365 break;
0366 }
0367 }
0368 }
0369 return phi_bin - 0.5;
0370 }
0371
0372 float PHG4TpcCylinderGeom::get_tbin_float(const double z) const
0373 {
0374 if (z < zmin || z > (zmin + nzbins * zstep))
0375 {
0376
0377 return -1;
0378 }
0379
0380 check_binning_method(PHG4CylinderCellDefs::sizebinning);
0381 return ((z - zmin) / zstep) - 0.5;
0382 }
0383
0384 int PHG4TpcCylinderGeom::get_phibin(const double phi, int side) const
0385 {
0386 double new_phi = phi;
0387 if (phi > M_PI)
0388 {
0389 new_phi = phi - 2 * M_PI;
0390 }
0391 if (phi < phimin)
0392 {
0393 new_phi = phi + 2 * M_PI;
0394 }
0395
0396 int phi_bin = find_phibin(new_phi, side);
0397
0398
0399
0400 if (phi_bin < 0)
0401 {
0402
0403 for (std::size_t s = 0; s < sector_max_Phi[side].size(); s++)
0404 {
0405 double daPhi = 0;
0406 if (s == 0)
0407 {
0408 daPhi = fabs(sector_min_Phi[side][11] + 2 * M_PI - sector_max_Phi[side][s]);
0409 }
0410 else
0411 {
0412 daPhi = fabs(sector_min_Phi[side][s - 1] - sector_max_Phi[side][s]);
0413 }
0414
0415 double min_phi = sector_max_Phi[side][s];
0416 double max_phi = sector_max_Phi[side][s] + daPhi;
0417 if (new_phi <= max_phi && new_phi >= min_phi)
0418 {
0419 if (fabs(max_phi - new_phi) > fabs(new_phi - min_phi))
0420 {
0421 new_phi = min_phi - phistep / 5;
0422 }
0423 else
0424 {
0425 new_phi = max_phi + phistep / 5;
0426 }
0427 }
0428 }
0429
0430
0431 phi_bin = find_phibin(new_phi, side);
0432 if (phi_bin < 0)
0433 {
0434 std::cout << PHWHERE << "Asking for bin for phi outside of phi range: " << phi << std::endl;
0435 exit(1);
0436
0437 }
0438 }
0439 return phi_bin;
0440 }
0441
0442 double
0443 PHG4TpcCylinderGeom::get_zcenter(const int ibin) const
0444 {
0445 if (ibin < 0 || ibin > nzbins)
0446 {
0447 std::cout << PHWHERE << "Asking for invalid bin in z: " << ibin << std::endl;
0448 exit(1);
0449 }
0450 check_binning_method(PHG4CylinderCellDefs::sizebinning);
0451 return zmin + (ibin + 0.5) * zstep;
0452 }
0453
0454 double
0455 PHG4TpcCylinderGeom::get_etacenter(const int ibin) const
0456 {
0457 if (ibin < 0 || ibin > nzbins)
0458 {
0459 std::cout << PHWHERE << "Asking for invalid bin in eta: " << ibin << std::endl;
0460 std::cout << "minbin: 0, maxbin " << nzbins << std::endl;
0461 exit(1);
0462 }
0463 check_binning_method_eta();
0464 return zmin + (ibin + 0.5) * zstep;
0465 }
0466
0467 double
0468 PHG4TpcCylinderGeom::get_phicenter_new(const int ibin) const
0469 {
0470 if (ibin < 0 || ibin > nphibins)
0471 {
0472 std::cout << PHWHERE << "Asking for invalid bin in phi: " << ibin << std::endl;
0473 exit(1);
0474 }
0475
0476 check_binning_method_phi();
0477
0478 return (phimin + (ibin + 0.5) * phistep);
0479 }
0480
0481 double
0482 PHG4TpcCylinderGeom::get_phicenter(const int ibin, const int side) const
0483 {
0484
0485 if (ibin < 0 || ibin > nphibins)
0486 {
0487 std::cout << PHWHERE << "Asking for invalid bin in phi: " << ibin << std::endl;
0488 exit(1);
0489 }
0490
0491 check_binning_method_phi();
0492
0493
0494 unsigned int pads_per_sector = nphibins / 12;
0495 unsigned int sector = ibin / pads_per_sector;
0496 double phi_center = (sector_max_Phi[side][sector] - (ibin + 0.5 - sector * pads_per_sector) * phistep);
0497 if (phi_center <= -M_PI)
0498 {
0499 phi_center += 2 * M_PI;
0500 }
0501 return phi_center;
0502 }
0503
0504 double
0505 PHG4TpcCylinderGeom::get_phi(const float ibin, const int side) const
0506 {
0507
0508 if (ibin < 0 || ibin > nphibins)
0509 {
0510 std::cout << PHWHERE << "Asking for invalid bin in phi: " << ibin << std::endl;
0511 exit(1);
0512 }
0513
0514 check_binning_method_phi();
0515
0516
0517 unsigned int pads_per_sector = nphibins / 12;
0518 unsigned int sector = ibin / pads_per_sector;
0519 double phi = (sector_max_Phi[side][sector] - (ibin + 0.5 - sector * pads_per_sector) * phistep);
0520 if (phi <= -M_PI)
0521 {
0522 phi += 2 * M_PI;
0523 }
0524 return phi;
0525 }
0526
0527 std::string
0528 PHG4TpcCylinderGeom::methodname(const int i) const
0529 {
0530 switch (i)
0531 {
0532 case PHG4CylinderCellDefs::sizebinning:
0533 return "Bins in cm";
0534 break;
0535 case PHG4CylinderCellDefs::etaphibinning:
0536 return "Eta/Phi bins";
0537 break;
0538 case PHG4CylinderCellDefs::etaslatbinning:
0539 return "Eta/numslat bins";
0540 break;
0541 case PHG4CylinderCellDefs::spacalbinning:
0542 return "SPACAL Tower bins";
0543 break;
0544 default:
0545 break;
0546 }
0547 return "Unknown";
0548 }
0549
0550 void PHG4TpcCylinderGeom::check_binning_method(const int i) const
0551 {
0552 if (binning != i)
0553 {
0554 std::cout << "different binning method used " << methodname(binning)
0555 << ", not : " << methodname(i)
0556 << std::endl;
0557 exit(1);
0558 }
0559 return;
0560 }
0561
0562 void PHG4TpcCylinderGeom::check_binning_method_eta(const std::string& src) const
0563 {
0564 if (binning != PHG4CylinderCellDefs::etaphibinning &&
0565 binning != PHG4CylinderCellDefs::etaslatbinning &&
0566 binning != PHG4CylinderCellDefs::spacalbinning)
0567 {
0568 if (src.size())
0569 {
0570 std::cout << src << " : ";
0571 }
0572
0573 std::cout << "different binning method used " << methodname(binning)
0574 << ", not : " << methodname(PHG4CylinderCellDefs::etaphibinning)
0575 << " or " << methodname(PHG4CylinderCellDefs::etaslatbinning)
0576 << " or " << methodname(PHG4CylinderCellDefs::spacalbinning)
0577 << std::endl;
0578 exit(1);
0579 }
0580 return;
0581 }
0582
0583 void PHG4TpcCylinderGeom::check_binning_method_phi(const std::string& src) const
0584 {
0585 if (binning != PHG4CylinderCellDefs::etaphibinning &&
0586 binning != PHG4CylinderCellDefs::sizebinning &&
0587 binning != PHG4CylinderCellDefs::etaslatbinning &&
0588 binning != PHG4CylinderCellDefs::spacalbinning)
0589 {
0590 if (src.size())
0591 {
0592 std::cout << src << " : ";
0593 }
0594
0595 std::cout << "different binning method used " << methodname(binning)
0596 << ", not : " << methodname(PHG4CylinderCellDefs::etaphibinning)
0597 << " or " << methodname(PHG4CylinderCellDefs::sizebinning)
0598 << " or " << methodname(PHG4CylinderCellDefs::etaslatbinning)
0599 << " or " << methodname(PHG4CylinderCellDefs::spacalbinning)
0600 << std::endl;
0601 exit(1);
0602 }
0603 return;
0604 }
