CaloGeometryTools/test/testCaloGeometryTools.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    testCaloCaloGeometryTools
0004 // Class:      testCaloCaloGeometryTools
0005 //
0006 /**\class testCaloCaloGeometryTools testCaloGeometryTools.cc test/testCaloGeometryTools/src/testCaloGeometryTools.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014
0015 // system include files
0016 #include <memory>
0017
0018 // user include files
0019 #include "FWCore/Framework/interface/stream/EDAnalyzer.h"
0020
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "FWCore/Framework/interface/Event.h"
0023 #include "FWCore/Framework/interface/ESHandle.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025
0026 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0027 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0028 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0029 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0030 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0031 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
0032 #include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
0033
0034 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0035 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0036
0037 #include "FastSimulation/CaloGeometryTools/interface/CaloGeometryHelper.h"
0038 #include "FastSimulation/CaloHitMakers/interface/EcalHitMaker.h"
0039 #include "FastSimulation/CaloGeometryTools/interface/Crystal.h"
0040 #include "FastSimulation/Utilities/interface/Histos.h"
0041 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0042
0043 #include <TH3F.h>
0044 #include <TPolyLine3D.h>
0045 #include <TMarker.h>
0046 #include <iostream>
0047 #include <sstream>
0048 //
0049 // class decleration
0050 //
0051
0052 typedef math::XYZVector XYZVector;
0053 typedef math::XYZVector XYZPoint;
0054
0055 class testCaloGeometryTools : public edm::stream::EDAnalyzer<> {
0056 public:
0057   explicit testCaloGeometryTools(const edm::ParameterSet&);
0058   ~testCaloGeometryTools();
0059
0060   virtual void analyze(const edm::Event&, const edm::EventSetup&);
0061
0062 private:
0063   // ----------member data ---------------------------
0064   void testpoint(const XYZPoint&, std::string name, bool barrel, RandomEngineAndDistribution const*);
0065   void checkSM();
0066   void checkSC();
0067   void testBorderCrossing();
0068
0069   const edm::ESGetToken<CaloTopology, CaloTopologyRecord> tokCaloTopo_;
0070   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tokCaloGeom_;
0071
0072   Histos* myHistos;
0073   CaloGeometryHelper myGeometry;
0074 };
0075
0076 //
0077 // constants, enums and typedefs
0078 //
0079
0080 //
0081 // static data member definitions
0082 //
0083
0084 //
0085 // constructors and destructor
0086 //
0087 testCaloGeometryTools::testCaloGeometryTools(const edm::ParameterSet& iConfig)
0088     : tokCaloTopo_(esConsumes()), tokCaloGeom_(esConsumes()) {
0089   myHistos = Histos::instance();
0090   myHistos->book("h100", 150, 0., 1.5, 100, 0., 35.);
0091 }
0092
0093 testCaloGeometryTools::~testCaloGeometryTools() {
0094   // do anything here that needs to be done at desctruction time
0095   // (e.g. close files, deallocate resources etc.)
0096   std::cout << " Writing the histo " << std::endl;
0097   myHistos->put("Grid.root");
0098   std::cout << " done " << std::endl;
0099 }
0100
0101 // ------------ method called to produce the data  ------------
0102 void testCaloGeometryTools::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0103   using namespace edm;
0104
0105   RandomEngineAndDistribution random(iEvent.streamID());
0106
0107   const edm::ESHandle<CaloTopology>& theCaloTopology = iSetup.getHandle(tokCaloTopo_);
0108   const edm::ESHandle<CaloGeometry>& pG = iSetup.getHandle(tokCaloGeom_);
0109
0110   // Setup the tools
0111   double bField000 = 4.;
0112   myGeometry.setupGeometry(*(pG.product()));
0113   myGeometry.setupTopology(*(theCaloTopology.product()));
0114   myGeometry.initialize(bField000);
0115
0116   // Take a point in the barrel
0117   XYZPoint p1(129, 0., -50);
0118   testpoint(p1, "barrel", true, &random);
0119   XYZPoint p2(60, 60, -317);
0120   testpoint(p1, "endcap", false, &random);
0121
0122   checkSM();
0123   checkSC();
0124   testBorderCrossing();
0125 }
0126
0127 void testCaloGeometryTools::checkSM() {
0128   const std::vector<DetId>& vec(myGeometry.getEcalBarrelGeometry()->getValidDetIds(DetId::Ecal, EcalBarrel));
0129   unsigned size = vec.size();
0130   for (unsigned ic = 0; ic < size; ++ic) {
0131     auto geom = (myGeometry.getEcalBarrelGeometry())->getGeometry(vec[ic]);
0132     GlobalPoint p = geom->getPosition();
0133     XYZPoint pp(p.x(), p.y(), p.z());
0134     // Build the name of the object
0135     std::ostringstream oss, oss2;
0136     oss << "iM" << ic;
0137     oss2 << "iSM" << ic;
0138     TMarker* myMarker = new TMarker(pp.eta(), pp.phi(), 1);
0139     myMarker->SetMarkerColor(EBDetId(vec[ic]).im());
0140     TMarker* myMarker2 = new TMarker(pp.eta(), pp.phi(), 1);
0141     myMarker2->SetMarkerColor(EBDetId(vec[ic]).ism());
0142     myHistos->addObject(oss.str(), myMarker);
0143     myHistos->addObject(oss2.str(), myMarker2);
0144   }
0145 }
0146
0147 void testCaloGeometryTools::checkSC() {
0148   const std::vector<DetId>& vec(myGeometry.getEcalEndcapGeometry()->getValidDetIds(DetId::Ecal, EcalEndcap));
0149   unsigned size = vec.size();
0150   for (unsigned ic = 0; ic < size; ++ic) {
0151     auto geom = (myGeometry.getEcalEndcapGeometry())->getGeometry(vec[ic]);
0152     GlobalPoint p = geom->getPosition();
0153     XYZPoint pp(p.x(), p.y(), p.z());
0154     // Build the name of the object
0155     std::ostringstream oss, oss2;
0156     if (p.z() > 0)
0157       oss << "iSCP" << ic;
0158     else
0159       oss << "iSCN" << ic;
0160     TMarker* myMarker = new TMarker(pp.x(), pp.y(), 1);
0161     if (pp.perp2() < 100.)
0162       std::cout << EEDetId(vec[ic]) << " " << pp.x() << " " << pp.y() << std::endl;
0163     myMarker->SetMarkerColor(EEDetId(vec[ic]).isc() % 100);
0164     myMarker->SetMarkerStyle(22);
0165     myHistos->addObject(oss.str(), myMarker);
0166   }
0167 }
0168
0169 void testCaloGeometryTools::testpoint(const XYZPoint& point,
0170                                       std::string name,
0171                                       bool barrel,
0172                                       RandomEngineAndDistribution const* random) {
0173   DetId myCell = myGeometry.getClosestCell(point, true, barrel);
0174   EcalHitMaker myGrid(&myGeometry, point, myCell, 1, 7, 0, random);
0175
0176   std::vector<Crystal> myCrystals = myGrid.getCrystals();
0177
0178   float xmin, ymin, zmin, xmax, ymax, zmax;
0179   xmin = ymin = zmin = 99999;
0180   xmax = ymax = zmax = -9999;
0181   unsigned nxtal = myCrystals.size();
0182
0183   std::vector<float> xp, yp, zp;
0184
0185   for (unsigned ic = 0; ic < nxtal; ++ic) {
0186     XYZPoint p = myCrystals[ic].getCenter();
0187
0188     myCrystals[ic].getDrawingCoordinates(xp, yp, zp);
0189     TPolyLine3D* myxtal = new TPolyLine3D(xp.size(), &xp[0], &yp[0], &zp[0]);
0190
0191     // Build the name of the object
0192     std::ostringstream oss;
0193     oss << name << ic;
0194
0195     myHistos->addObject(oss.str(), myxtal);
0196
0197     if (xmin > p.x())
0198       xmin = p.x();
0199     if (ymin > p.y())
0200       ymin = p.y();
0201     if (zmin > p.z())
0202       zmin = p.z();
0203     if (xmax < p.x())
0204       xmax = p.x();
0205     if (ymax < p.y())
0206       ymax = p.y();
0207     if (zmax < p.z())
0208       zmax = p.z();
0209   }
0210   TH3F* frame = new TH3F(std::string(name + "frame").c_str(),
0211                          "",
0212                          100,
0213                          xmin * 0.9,
0214                          xmax * 1.1,
0215                          100,
0216                          ymin * 0.9,
0217                          ymax * 1.1,
0218                          100,
0219                          zmin * 0.9,
0220                          zmax * 1.1);
0221   myHistos->addObject("frame" + name, frame);
0222 }
0223
0224 void testCaloGeometryTools::testBorderCrossing() {
0225   // Barrel
0226   const std::vector<DetId>& vec(myGeometry.getEcalBarrelGeometry()->getValidDetIds(DetId::Ecal, EcalBarrel));
0227   unsigned size = vec.size();
0228   unsigned counter = 0;
0229   for (unsigned ic = 0; ic < size; ++ic) {
0230     CaloGeometryHelper::NeiVect neighbours = myGeometry.getNeighbours(vec[ic]);
0231     for (unsigned in = 0; in < 8; ++in) {
0232       if (neighbours[in].null())
0233         continue;
0234       if (myGeometry.borderCrossing(vec[ic], neighbours[in])) {
0235         auto geom = (myGeometry.getEcalBarrelGeometry())->getGeometry(vec[ic]);
0236         GlobalPoint p1 = geom->getPosition();
0237         XYZPoint pp1(p1.x(), p1.y(), p1.z());
0238         geom = (myGeometry.getEcalBarrelGeometry())->getGeometry(neighbours[in]);
0239         GlobalPoint p2 = geom->getPosition();
0240         XYZPoint pp2(p2.x(), p2.y(), p2.z());
0241         TMarker* myMarker = new TMarker((pp1 + pp2).eta() * 0.5, ((pp1 + pp2) * 0.5).phi(), 22);
0242         std::ostringstream oss;
0243         oss << "iBCB" << counter;
0244         myHistos->addObject(oss.str(), myMarker);
0245         ++counter;
0246       }
0247     }
0248   }
0249
0250   // Endcap
0251   const std::vector<DetId>& vec2(myGeometry.getEcalEndcapGeometry()->getValidDetIds(DetId::Ecal, EcalEndcap));
0252   size = vec2.size();
0253   counter = 0;
0254   for (unsigned ic = 0; ic < size; ++ic) {
0255     CaloGeometryHelper::NeiVect neighbours = myGeometry.getNeighbours(vec2[ic]);
0256     for (unsigned in = 0; in < 8; ++in) {
0257       if (neighbours[in].null())
0258         continue;
0259       if (myGeometry.borderCrossing(vec2[ic], neighbours[in])) {
0260         auto geom = (myGeometry.getEcalEndcapGeometry())->getGeometry(vec2[ic]);
0261         GlobalPoint p1 = geom->getPosition();
0262         XYZPoint pp1(p1.x(), p1.y(), p1.z());
0263         geom = (myGeometry.getEcalEndcapGeometry())->getGeometry(neighbours[in]);
0264         GlobalPoint p2 = geom->getPosition();
0265         XYZPoint pp2(p2.x(), p2.y(), p2.z());
0266         TMarker* myMarker = new TMarker((pp1 + pp2).x() * 0.5, (pp1 + pp2).y() * 0.5, 22);
0267         std::ostringstream oss;
0268         if (p1.z() > 0)
0269           oss << "iBCEN" << counter;
0270         else
0271           oss << "iBCEP" << counter;
0272
0273         myHistos->addObject(oss.str(), myMarker);
0274         ++counter;
0275       }
0276     }
0277   }
0278 }
0279
0280 //define this as a plug-in
0281
0282 DEFINE_FWK_MODULE(testCaloGeometryTools);