Back to home page

Project CMSSW displayed by LXR
 
 

    

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) {  // don't copy into yourself
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   //  std::cout << " Hello " << std::endl;
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     // Set explicity the z to 0 !
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   //  std::cout << " End of 1 constructor " << std::endl;
0046   //  std::cout << " Ncorners " << corners_.size() << std::endl;
0047   //  std::cout << " Ndirs " << dir_.size() << std::endl;
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   //  std::cout << " We are in the 2nd constructor " << std::endl;
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     // the good one in the central
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     //      std::cout << " Constructor 2; input corners "  << std::endl;
0074     for (unsigned ic = 0; ic < 4; ++ic) {
0075       //      std::cout << corners[ic]<< " " ;
0076       XYZPoint corner = rotation_(corners[ic]) + translation_;
0077       //      std::cout << corner << std::endl ;
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   //  std::cout << " End of 2 constructor " << std::endl;
0087   //  std::cout << " Corners(constructor) " ;
0088   //  std::cout << corners_[0] << std::endl;
0089   //  std::cout << corners_[1] << std::endl;
0090   //  std::cout << corners_[2] << std::endl;
0091   //  std::cout << corners_[3] << std::endl;
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   //  std::cout << " We are in the 2nd constructor " << std::endl;
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     // the good one in the central
0110     trans_ = trans;
0111     //      std::cout << " Constructor 2; input corners "  << std::endl;
0112     trans_.GetDecomposition(rotation_, translation_);
0113     for (unsigned ic = 0; ic < 4; ++ic) {
0114       XYZPoint corner = rotation_(corners[ic]) + translation_;
0115       //      std::cout << corner << std::endl ;
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   //  std::cout << "Inside " << ppoint <<std::endl;
0129   //  std::cout << "Corners " << corners_.size() << std::endl;
0130   //  std::cout << corners_[0] << std::endl;
0131   //  std::cout << corners_[1] << std::endl;
0132   //  std::cout << corners_[2] << std::endl;
0133   //  std::cout << corners_[3] << std::endl;
0134   //  std::cout << " Got the 2D point " << std::endl;
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   //  double N1(n1.mag());
0141   //  double N2(n2.mag());
0142   double r1(n1 * n2);
0143   bool inside1(r1 <= 0.);
0144 
0145   if (!inside1)
0146     return false;
0147 
0148   //  if(debug)
0149   //    {
0150   //      std::cout << n1 << std::endl;
0151   //      std::cout << n2 << std::endl;
0152   //      std::cout << r1 << std::endl;
0153   //      std::cout << inside1 << std::endl;
0154   //    }
0155 
0156   //  bool close1=(N1<epsilon_||N2<epsilon_);
0157   //
0158   //  if(!close1&&!inside1) return false;
0159   //  std::cout << " First calculation " << std::endl;
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   //  double N3(n3.mag());
0165   //  double N4(n4.mag());
0166   //  bool close2=(N3<epsilon_||N4<epsilon_);
0167   double r2(n3 * n4);
0168   bool inside2(r2 <= 0.);
0169   //  //  std::cout << " pv1 & pv3 " << pv1.mag() << " " << pv3.mag() << std::endl;
0170   //  //  double tmp=(pv1-(pv1*dir_[1])*dir_[1])*(pv3-(pv3*dir_[3])*dir_[3]);
0171   //  //  std::cout << " Computed tmp " << tmp << std::endl;
0172   //  if(debug)
0173   //    {
0174   //      std::cout << n3 << std::endl;
0175   //      std::cout << n4 << std::endl;
0176   //      std::cout << r2 << std::endl;
0177   //      std::cout << inside2 << std::endl;
0178   //    }
0179   //  if(!close2&&!inside2) return false;
0180   //  std::cout << " Second calculation " << std::endl;
0181   //  std::cout << " True " << std::endl;
0182   //    return (!close1&&!close2||(close2&&inside1||close1&&inside2));
0183 
0184   return inside2;
0185 }
0186 
0187 /*
0188 bool 
0189 CrystalPad::globalinside(XYZPoint point) const
0190 {
0191   //  std::cout << " Global inside " << std::endl;
0192   //  std::cout << point << " " ;
0193   ROOT::Math::Rotation3D r;
0194   XYZVector t;
0195   point = rotation_(point)+translation_;
0196   //  std::cout << point << std::endl;
0197   //  print();
0198   CLHEP::Hep2Vector ppoint(point.X(),point.Y());
0199   bool result=inside(ppoint);
0200   //  std::cout << " Result " << result << std::endl;
0201   return result;
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 CLHEP::Hep2Vector 
0215 CrystalPad::localPoint(XYZPoint point) const
0216 {
0217   point = rotation_(point)+translation_;
0218   return CLHEP::Hep2Vector(point.X(),point.Y());
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   // Find the centre-of-gravity of the Quad (after re-organization)
0263   center_ = CLHEP::Hep2Vector(0., 0.);
0264   for (unsigned ic = 0; ic < 4; ++ic)
0265     center_ += corners_[ic];
0266   center_ *= 0.25;
0267 
0268   // Rescale the corners to allow for some inaccuracies in
0269   // in the inside test
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!