EcalTestBeam/test/CrystalCenterDump.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    CrystalCenterDump
0004 // Class:      CrystalCenterDump
0005 //
0006 /**\class CrystalCenterDump CrystalCenterDump.cc test/CrystalCenterDump/src/CrystalCenterDump.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 #include <cmath>
0018 #include <iomanip>
0019
0020 // user include files
0021 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0022
0023 #include "FWCore/Framework/interface/EventSetup.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0026 #include "FWCore/Utilities/interface/ESGetToken.h"
0027
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0030 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0031 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0032 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0033 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0034 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0035 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0036
0037 #include <fstream>
0038
0039 #include "CLHEP/Vector/ThreeVector.h"
0040
0041 class CrystalCenterDump : public edm::one::EDAnalyzer<> {
0042 public:
0043   explicit CrystalCenterDump(const edm::ParameterSet&);
0044   ~CrystalCenterDump() override;
0045
0046   void beginJob() override {}
0047   void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0048   void endJob() override {}
0049
0050 private:
0051   void build(const CaloGeometry& cg, DetId::Detector det, int subdetn, const char* name);
0052   int pass_;
0053
0054   double A_;
0055   double B_;
0056   double beamEnergy_;
0057
0058   double crystalDepth() { return A_ * (B_ + log(beamEnergy_)); }
0059
0060   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geometryToken_;
0061 };
0062
0063 CrystalCenterDump::CrystalCenterDump(const edm::ParameterSet& iConfig) {
0064   //now do what ever initialization is needed
0065   pass_ = 0;
0066
0067   A_ = iConfig.getUntrackedParameter<double>("Afac", 0.89);
0068   B_ = iConfig.getUntrackedParameter<double>("Bfac", 5.7);
0069   beamEnergy_ = iConfig.getUntrackedParameter<double>("BeamEnergy", 120.);
0070
0071   edm::LogVerbatim("CrysInfo") << "Position computed according to the depth " << crystalDepth() << " based on:"
0072                                << "\n A = " << A_ << " cm "
0073                                << "\n B = " << B_ << "\n BeamEnergy = " << beamEnergy_ << " GeV";
0074
0075   geometryToken_ = esConsumes<CaloGeometry, CaloGeometryRecord>(edm::ESInputTag{});
0076 }
0077
0078 CrystalCenterDump::~CrystalCenterDump() {
0079   // do anything here that needs to be done at desctruction time
0080   // (e.g. close files, deallocate resources etc.)
0081 }
0082
0083 void CrystalCenterDump::build(const CaloGeometry& cg, DetId::Detector det, int subdetn, const char* name) {
0084   std::fstream f(name, std::ios_base::out);
0085   const CaloSubdetectorGeometry* geom(cg.getSubdetectorGeometry(det, subdetn));
0086
0087   const std::vector<DetId>& ids = geom->getValidDetIds(det, subdetn);
0088   for (auto id : ids) {
0089     auto cell = geom->getGeometry(id);
0090     if (det == DetId::Ecal) {
0091       if (subdetn == EcalBarrel) {
0092         EBDetId ebid(id.rawId());
0093         if (ebid.ism() == 1) {
0094           float depth = (crystalDepth());
0095           double crysX = (cell)->getPosition(depth).x();
0096           double crysY = (cell)->getPosition(depth).y();
0097           double crysZ = (cell)->getPosition(depth).z();
0098
0099           CLHEP::Hep3Vector crysPos(crysX, crysY, crysZ);
0100           double crysEta = crysPos.eta();
0101           double crysTheta = crysPos.theta();
0102           double crysPhi = crysPos.phi();
0103
0104           edm::LogVerbatim("CrysPos") << ebid.ic() << " x = " << crysX << " y = " << crysY << " z = " << crysZ << " \n "
0105                                       << " eta = " << crysEta << " phi = " << crysPhi << " theta = " << crysTheta;
0106           f << std::setw(4) << ebid.ic() << " " << std::setw(8) << std::setprecision(6) << crysEta << " "
0107             << std::setw(8) << std::setprecision(6) << crysPhi << " " << std::endl;
0108         }
0109       }
0110     }
0111   }
0112   f.close();
0113 }
0114
0115 //
0116 // member functions
0117 //
0118
0119 // ------------ method called to produce the data  ------------
0120 void CrystalCenterDump::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0121   edm::LogVerbatim("CrysPos") << "Writing the center (eta,phi) for crystals in barrel SM 1 ";
0122
0123   const auto& pG = iSetup.getData(geometryToken_);
0124   //
0125   // get the ecal & hcal geometry
0126   //
0127   if (pass_ == 0) {
0128     build(pG, DetId::Ecal, EcalBarrel, "BarrelSM1CrystalCenter.dat");
0129   }
0130
0131   pass_++;
0132 }
0133
0134 //define this as a plug-in
0135
0136 DEFINE_FWK_MODULE(CrystalCenterDump);