File indexing completed on 2021-02-14 13:08:05
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 #include <memory>
0017 #include <cmath>
0018 #include <iomanip>
0019
0020
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
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
0080
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
0119
0120
0121
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
0128
0129 if (pass_ == 0) {
0130 build(pG, DetId::Ecal, EcalBarrel, "BarrelSM1CrystalCenter.dat");
0131 }
0132
0133 pass_++;
0134 }
0135
0136
0137
0138 DEFINE_FWK_MODULE(CrystalCenterDump);