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::LogInfo("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   int n = 0;
0088   const std::vector<DetId>& ids = geom->getValidDetIds(det, subdetn);
0089   for (auto id : ids) {
0090     n++;
0091     auto cell = geom->getGeometry(id);
0092     if (det == DetId::Ecal) {
0093       if (subdetn == EcalBarrel) {
0094         EBDetId ebid(id.rawId());
0095         if (ebid.ism() == 1) {
0096           float depth = (crystalDepth());
0097           double crysX = (cell)->getPosition(depth).x();
0098           double crysY = (cell)->getPosition(depth).y();
0099           double crysZ = (cell)->getPosition(depth).z();
0100
0101           CLHEP::Hep3Vector crysPos(crysX, crysY, crysZ);
0102           double crysEta = crysPos.eta();
0103           double crysTheta = crysPos.theta();
0104           double crysPhi = crysPos.phi();
0105
0106           edm::LogInfo("CrysPos") << ebid.ic() << " x = " << crysX << " y = " << crysY << " z = " << crysZ << " \n "
0107                                   << " eta = " << crysEta << " phi = " << crysPhi << " theta = " << crysTheta;
0108           f << std::setw(4) << ebid.ic() << " " << std::setw(8) << std::setprecision(6) << crysEta << " "
0109             << std::setw(8) << std::setprecision(6) << crysPhi << " " << std::endl;
0110         }
0111       }
0112     }
0113   }
0114   f.close();
0115 }
0116
0117 //
0118 // member functions
0119 //
0120
0121 // ------------ method called to produce the data  ------------
0122 void CrystalCenterDump::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0123   std::cout << "Writing the center (eta,phi) for crystals in barrel SM 1 " << std::endl;
0124
0125   const auto& pG = iSetup.getData(geometryToken_);
0126   //
0127   // get the ecal & hcal geometry
0128   //
0129   if (pass_ == 0) {
0130     build(pG, DetId::Ecal, EcalBarrel, "BarrelSM1CrystalCenter.dat");
0131   }
0132
0133   pass_++;
0134 }
0135
0136 //define this as a plug-in
0137
0138 DEFINE_FWK_MODULE(CrystalCenterDump);