File indexing completed on 2023-10-25 09:45:59
0001 #include "FastSimulation/CaloGeometryTools/interface/CrystalPad.h"
0002
0003 #include <iostream>
0004
0005 std::vector<CLHEP::Hep2Vector> CrystalPad::aVector(4);
0006
0007 CrystalPad::CrystalPad(const CrystalPad& right) {
0008 corners_ = right.corners_;
0009 dir_ = right.dir_;
0010 number_ = right.number_;
0011 survivalProbability_ = right.survivalProbability_;
0012 center_ = right.center_;
0013 epsilon_ = right.epsilon_;
0014 dummy_ = right.dummy_;
0015 }
0016
0017 CrystalPad& CrystalPad::operator=(const CrystalPad& right) {
0018 if (this != &right) {
0019 corners_ = right.corners_;
0020 dir_ = right.dir_;
0021 number_ = right.number_;
0022 survivalProbability_ = right.survivalProbability_;
0023 center_ = right.center_;
0024 epsilon_ = right.epsilon_;
0025 dummy_ = right.dummy_;
0026 }
0027 return *this;
0028 }
0029
0030 CrystalPad::CrystalPad(unsigned number, const std::vector<CLHEP::Hep2Vector>& corners)
0031 : corners_(corners), dir_(aVector), number_(number), survivalProbability_(1.), center_(0., 0.), epsilon_(0.001) {
0032
0033 if (corners.size() != 4) {
0034 std::cout << " Try to construct a quadrilateral with " << corners.size() << " points ! " << std::endl;
0035 dummy_ = true;
0036 } else {
0037 dummy_ = false;
0038
0039 for (unsigned ic = 0; ic < 4; ++ic) {
0040 dir_[ic] = (corners[(ic + 1) % 4] - corners[ic]).unit();
0041 center_ += corners_[ic];
0042 }
0043 center_ *= 0.25;
0044 }
0045
0046
0047
0048 }
0049
0050 CrystalPad::CrystalPad(unsigned number,
0051 int onEcal,
0052 const std::vector<XYZPoint>& corners,
0053 const XYZPoint& origin,
0054 const XYZVector& vec1,
0055 const XYZVector& vec2)
0056 : corners_(aVector), dir_(aVector), number_(number), survivalProbability_(1.), center_(0., 0.), epsilon_(0.001) {
0057
0058 if (corners.size() != 4) {
0059 std::cout << " Try to construct a quadrilateral with " << corners.size() << " points ! " << std::endl;
0060 dummy_ = true;
0061 } else {
0062 dummy_ = false;
0063 double sign = (onEcal == 1) ? -1. : 1.;
0064
0065
0066 trans_ = Transform3D((Point)origin,
0067 (Point)(origin + vec1),
0068 (Point)(origin + vec2),
0069 Point(0., 0., 0.),
0070 Point(0., 0., sign),
0071 Point(0., 1., 0.));
0072 trans_.GetDecomposition(rotation_, translation_);
0073
0074 for (unsigned ic = 0; ic < 4; ++ic) {
0075
0076 XYZPoint corner = rotation_(corners[ic]) + translation_;
0077
0078 corners_[ic] = CLHEP::Hep2Vector(corner.X(), corner.Y());
0079 center_ += corners_[ic];
0080 }
0081 for (unsigned ic = 0; ic < 4; ++ic) {
0082 dir_[ic] = (corners_[(ic + 1) % 4] - corners_[ic]).unit();
0083 }
0084 center_ *= 0.25;
0085 }
0086
0087
0088
0089
0090
0091
0092 }
0093 CrystalPad::CrystalPad(
0094 unsigned number, const std::vector<XYZPoint>& corners, const Transform3D& trans, double scaf, bool bothdirections)
0095 : corners_(aVector),
0096 dir_(aVector),
0097 number_(number),
0098 survivalProbability_(1.),
0099 center_(0., 0.),
0100 epsilon_(0.001),
0101 yscalefactor_(scaf) {
0102
0103 if (corners.size() != 4) {
0104 std::cout << " Try to construct a quadrilateral with " << corners.size() << " points ! " << std::endl;
0105 dummy_ = true;
0106 } else {
0107 dummy_ = false;
0108
0109
0110 trans_ = trans;
0111
0112 trans_.GetDecomposition(rotation_, translation_);
0113 for (unsigned ic = 0; ic < 4; ++ic) {
0114 XYZPoint corner = rotation_(corners[ic]) + translation_;
0115
0116 double xscalefactor = (bothdirections) ? yscalefactor_ : 1.;
0117 corners_[ic] = CLHEP::Hep2Vector(corner.X() * xscalefactor, corner.Y() * yscalefactor_);
0118 center_ += corners_[ic];
0119 }
0120 for (unsigned ic = 0; ic < 4; ++ic) {
0121 dir_[ic] = (corners_[(ic + 1) % 4] - corners_[ic]).unit();
0122 }
0123 center_ *= 0.25;
0124 }
0125 }
0126
0127 bool CrystalPad::inside(const CLHEP::Hep2Vector& ppoint, bool debug) const {
0128
0129
0130
0131
0132
0133
0134
0135 CLHEP::Hep2Vector pv0(ppoint - corners_[0]);
0136 CLHEP::Hep2Vector pv2(ppoint - corners_[2]);
0137 CLHEP::Hep2Vector n1(pv0 - (pv0 * dir_[0]) * dir_[0]);
0138 CLHEP::Hep2Vector n2(pv2 - (pv2 * dir_[2]) * dir_[2]);
0139
0140
0141
0142 double r1(n1 * n2);
0143 bool inside1(r1 <= 0.);
0144
0145 if (!inside1)
0146 return false;
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160 CLHEP::Hep2Vector pv1(ppoint - corners_[1]);
0161 CLHEP::Hep2Vector pv3(ppoint - corners_[3]);
0162 CLHEP::Hep2Vector n3(pv1 - (pv1 * dir_[1]) * dir_[1]);
0163 CLHEP::Hep2Vector n4(pv3 - (pv3 * dir_[3]) * dir_[3]);
0164
0165
0166
0167 double r2(n3 * n4);
0168 bool inside2(r2 <= 0.);
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184 return inside2;
0185 }
0186
0187
0188
0189
0190
0191
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201
0202
0203
0204
0205 void CrystalPad::print() const {
0206 std::cout << " Corners " << std::endl;
0207 std::cout << corners_[0] << std::endl;
0208 std::cout << corners_[1] << std::endl;
0209 std::cout << corners_[2] << std::endl;
0210 std::cout << corners_[3] << std::endl;
0211 }
0212
0213
0214
0215
0216
0217
0218
0219
0220
0221
0222 CLHEP::Hep2Vector& CrystalPad::edge(unsigned iside, int n) { return corners_[(iside + n) % 4]; }
0223
0224 CLHEP::Hep2Vector& CrystalPad::edge(CaloDirection dir) {
0225 switch (dir) {
0226 case NORTHWEST:
0227 return corners_[0];
0228 break;
0229 case NORTHEAST:
0230 return corners_[1];
0231 break;
0232 case SOUTHEAST:
0233 return corners_[2];
0234 break;
0235 case SOUTHWEST:
0236 return corners_[3];
0237 break;
0238 default: {
0239 std::cout << " Serious problem in CrystalPad ! " << dir << std::endl;
0240 return corners_[0];
0241 }
0242 }
0243 return corners_[0];
0244 }
0245
0246 void CrystalPad::extrems(double& xmin, double& xmax, double& ymin, double& ymax) const {
0247 xmin = ymin = 999;
0248 xmax = ymax = -999;
0249 for (unsigned ic = 0; ic < 4; ++ic) {
0250 if (corners_[ic].x() < xmin)
0251 xmin = corners_[ic].x();
0252 if (corners_[ic].x() > xmax)
0253 xmax = corners_[ic].x();
0254 if (corners_[ic].y() < ymin)
0255 ymin = corners_[ic].y();
0256 if (corners_[ic].y() > ymax)
0257 ymax = corners_[ic].y();
0258 }
0259 }
0260
0261 void CrystalPad::resetCorners() {
0262
0263 center_ = CLHEP::Hep2Vector(0., 0.);
0264 for (unsigned ic = 0; ic < 4; ++ic)
0265 center_ += corners_[ic];
0266 center_ *= 0.25;
0267
0268
0269
0270 for (unsigned ic = 0; ic < 4; ++ic)
0271 corners_[ic] += 0.001 * (corners_[ic] - center_);
0272 }
0273
0274 std::ostream& operator<<(std::ostream& ost, CrystalPad& quad) {
0275 ost << " Number " << quad.getNumber() << std::endl;
0276 ost << NORTHWEST << quad.edge(NORTHWEST) << std::endl;
0277 ost << NORTHEAST << quad.edge(NORTHEAST) << std::endl;
0278 ost << SOUTHEAST << quad.edge(SOUTHEAST) << std::endl;
0279 ost << SOUTHWEST << quad.edge(SOUTHWEST) << std::endl;
0280
0281 return ost;
0282 }
0283
0284 void CrystalPad::getDrawingCoordinates(std::vector<float>& x, std::vector<float>& y) const {
0285 x.clear();
0286 y.clear();
0287 x.push_back(corners_[0].x());
0288 x.push_back(corners_[1].x());
0289 x.push_back(corners_[2].x());
0290 x.push_back(corners_[3].x());
0291 x.push_back(corners_[0].x());
0292 y.push_back(corners_[0].y());
0293 y.push_back(corners_[1].y());
0294 y.push_back(corners_[2].y());
0295 y.push_back(corners_[3].y());
0296 y.push_back(corners_[0].y());
0297 }