File indexing completed on 2024-09-11 04:33:23
0001
0002
0003
0004
0005
0006 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0007 #include "SimG4Core/Watcher/interface/SimProducer.h"
0008 #include "SimG4Core/Notification/interface/Observer.h"
0009
0010
0011 #include "DataFormats/DetId/interface/DetId.h"
0012 #include "DataFormats/Math/interface/Point3D.h"
0013 #include "SimDataFormats/ValidationFormats/interface/PHGCalValidInfo.h"
0014 #include "SimG4CMS/Calo/interface/HGCalNumberingScheme.h"
0015
0016 #include "FWCore/Framework/interface/Event.h"
0017 #include "FWCore/Framework/interface/EventSetup.h"
0018 #include "FWCore/Framework/interface/ESHandle.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021
0022 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0023 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0024
0025 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0026 #include "FWCore/PluginManager/interface/ModuleDef.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028
0029 #include "G4Step.hh"
0030 #include "G4Track.hh"
0031 #include "G4NavigationHistory.hh"
0032 #include "G4TouchableHistory.hh"
0033
0034 #include <CLHEP/Units/SystemOfUnits.h>
0035 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0036
0037 #include <cmath>
0038 #include <iostream>
0039 #include <iomanip>
0040 #include <memory>
0041 #include <vector>
0042 #include <string>
0043
0044 class SimG4HGCalValidation : public SimProducer, public Observer<const BeginOfEvent*>, public Observer<const G4Step*> {
0045 public:
0046 SimG4HGCalValidation(const edm::ParameterSet& p);
0047 ~SimG4HGCalValidation() override = default;
0048
0049 void registerConsumes(edm::ConsumesCollector) override;
0050 void produce(edm::Event&, const edm::EventSetup&) override;
0051 void beginRun(edm::EventSetup const&) override;
0052
0053 private:
0054 SimG4HGCalValidation(const SimG4HGCalValidation&);
0055 const SimG4HGCalValidation& operator=(const SimG4HGCalValidation&);
0056
0057 void init();
0058
0059
0060 void update(const BeginOfEvent* evt) override;
0061 void update(const G4Step* step) override;
0062
0063
0064 void layerAnalysis(PHGCalValidInfo&);
0065 void clear();
0066
0067 private:
0068
0069 std::vector<edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord>> ddconsToken_;
0070 std::vector<std::unique_ptr<HGCalNumberingScheme>> hgcalNumbering_;
0071
0072
0073 const edm::ParameterSet m_Anal;
0074 const std::vector<std::string> names_;
0075 const std::vector<int> types_, detTypes_;
0076 const std::string labelLayer_;
0077 const int verbosity_;
0078 std::vector<int> subdet_;
0079 std::vector<std::string> nameXs_;
0080
0081
0082 int levelT1_, levelT2_;
0083
0084
0085 unsigned int count_;
0086 double edepEE_, edepHEF_, edepHEB_;
0087 std::vector<double> hgcEEedep_, hgcHEFedep_, hgcHEBedep_;
0088 std::vector<unsigned int> dets_, hgchitDets_, hgchitIndex_;
0089 std::vector<double> hgchitX_, hgchitY_, hgchitZ_;
0090 };
0091
0092 SimG4HGCalValidation::SimG4HGCalValidation(const edm::ParameterSet& p)
0093 : m_Anal(p.getParameter<edm::ParameterSet>("SimG4HGCalValidation")),
0094 names_(m_Anal.getParameter<std::vector<std::string>>("Names")),
0095 types_(m_Anal.getParameter<std::vector<int>>("Types")),
0096 detTypes_(m_Anal.getParameter<std::vector<int>>("DetTypes")),
0097 labelLayer_(m_Anal.getParameter<std::string>("LabelLayerInfo")),
0098 verbosity_(m_Anal.getUntrackedParameter<int>("Verbosity", 0)),
0099 levelT1_(999),
0100 levelT2_(999),
0101 count_(0) {
0102 produces<PHGCalValidInfo>(labelLayer_);
0103
0104 if (verbosity_ > 0) {
0105 edm::LogVerbatim("ValidHGCal") << "HGCalTestAnalysis:: Initialised as "
0106 << "observer of begin events and of G4step "
0107 << "with Parameter values: \n\tLabel : " << labelLayer_ << " and with "
0108 << names_.size() << " detectors";
0109 for (unsigned int k = 0; k < names_.size(); ++k)
0110 edm::LogVerbatim("ValidHGCal") << " [" << k << "] " << names_[k] << " Type " << types_[k] << " DetType "
0111 << detTypes_[k];
0112 }
0113 }
0114
0115 void SimG4HGCalValidation::registerConsumes(edm::ConsumesCollector cc) {
0116 for (unsigned int type = 0; type < types_.size(); ++type) {
0117 int detType = detTypes_[type];
0118 G4String nameX = "HGCal";
0119 if (types_[type] <= 1) {
0120 dets_.emplace_back(static_cast<unsigned int>(DetId::Forward));
0121 if (detType == 0)
0122 subdet_.emplace_back(static_cast<int>(ForwardSubdetector::HGCEE));
0123 else if (detType == 1)
0124 subdet_.emplace_back(static_cast<int>(ForwardSubdetector::HGCHEF));
0125 else
0126 subdet_.emplace_back(static_cast<int>(ForwardSubdetector::HGCHEB));
0127 if (detType == 0)
0128 nameX = "HGCalEESensitive";
0129 else if (detType == 1)
0130 nameX = "HGCalHESiliconSensitive";
0131 else
0132 nameX = "HGCalHEScintillatorSensitive";
0133 nameXs_.emplace_back(nameX);
0134 ddconsToken_.emplace_back(
0135 cc.esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag{"", nameX}));
0136 edm::LogVerbatim("ValidHGCal") << "Defines ESGetToken for type " << type << ":" << dets_.back() << ":"
0137 << subdet_.back() << ":" << nameX;
0138 } else {
0139 edm::LogError("ValidHGCal") << "Wrong Type " << types_[type];
0140 throw cms::Exception("Unknown", "ValidHGCal") << "Wrong Type " << types_[type] << "\n";
0141 }
0142 }
0143 }
0144
0145 void SimG4HGCalValidation::produce(edm::Event& e, const edm::EventSetup&) {
0146 std::unique_ptr<PHGCalValidInfo> productLayer(new PHGCalValidInfo);
0147 layerAnalysis(*productLayer);
0148 e.put(std::move(productLayer), labelLayer_);
0149 }
0150
0151 void SimG4HGCalValidation::beginRun(edm::EventSetup const& es) {
0152 for (unsigned int type = 0; type < types_.size(); ++type) {
0153 int layers(0);
0154 int detType = detTypes_[type];
0155 const edm::ESHandle<HGCalDDDConstants>& hdc = es.getHandle(ddconsToken_[type]);
0156 if (hdc.isValid()) {
0157 levelT1_ = hdc->levelTop(0);
0158 levelT2_ = hdc->levelTop(1);
0159 hgcalNumbering_.emplace_back(
0160 std::make_unique<HGCalNumberingScheme>(*hdc, static_cast<DetId::Detector>(dets_[type]), nameXs_[type], ""));
0161 layers = hdc->layers(false);
0162 } else {
0163 edm::LogError("ValidHGCal") << "Cannot find HGCalDDDConstants for " << nameXs_[type];
0164 throw cms::Exception("Unknown", "ValidHGCal") << "Cannot find HGCalDDDConstants for " << nameXs_[type] << "\n";
0165 }
0166 if (detType == 0) {
0167 for (int i = 0; i < layers; ++i)
0168 hgcEEedep_.emplace_back(0);
0169 } else if (detType == 1) {
0170 for (int i = 0; i < layers; ++i)
0171 hgcHEFedep_.emplace_back(0);
0172 } else {
0173 for (int i = 0; i < layers; ++i)
0174 hgcHEBedep_.emplace_back(0);
0175 }
0176 if (verbosity_ > 0)
0177 edm::LogVerbatim("ValidHGCal") << "[" << type << "]: " << nameXs_[type] << " det " << dets_[type] << " subdet "
0178 << subdet_[type] << " with " << layers << " layers";
0179 }
0180 }
0181
0182
0183 void SimG4HGCalValidation::update(const BeginOfEvent* evt) {
0184 int iev = (*evt)()->GetEventID();
0185 if (verbosity_ > 0)
0186 edm::LogVerbatim("ValidHGCal") << "SimG4HGCalValidation: =====> Begin "
0187 << "event = " << iev;
0188
0189 ++count_;
0190 edepEE_ = edepHEF_ = edepHEB_ = 0.;
0191
0192
0193 for (unsigned int i = 0; i < hgcEEedep_.size(); i++)
0194 hgcEEedep_[i] = 0.;
0195 for (unsigned int i = 0; i < hgcHEFedep_.size(); i++)
0196 hgcHEFedep_[i] = 0.;
0197 for (unsigned int i = 0; i < hgcHEBedep_.size(); i++)
0198 hgcHEBedep_[i] = 0.;
0199
0200
0201 clear();
0202 }
0203
0204
0205 void SimG4HGCalValidation::update(const G4Step* aStep) {
0206 if (aStep != nullptr) {
0207 G4VPhysicalVolume* curPV = aStep->GetPreStepPoint()->GetPhysicalVolume();
0208 const G4String& name = curPV->GetName();
0209 G4VSensitiveDetector* curSD = aStep->GetPreStepPoint()->GetSensitiveDetector();
0210 if (verbosity_ > 1)
0211 edm::LogVerbatim("ValidHGCal") << "ValidHGCal::Step in " << name << " at "
0212 << aStep->GetPreStepPoint()->GetPosition();
0213
0214
0215 if (curSD != nullptr) {
0216 int type(-1);
0217 for (unsigned int k = 0; k < names_.size(); ++k) {
0218 if (name.find(names_[k].c_str()) != std::string::npos) {
0219 type = k;
0220 break;
0221 }
0222 }
0223 int detType = (type >= 0) ? detTypes_[type] : -1;
0224 if (verbosity_ > 0)
0225 edm::LogVerbatim("ValidHGCal") << "ValidHGCal::In SD " << name << " type " << type << ":" << detType;
0226
0227
0228 if (type >= 0) {
0229
0230 const G4TouchableHistory* touchable =
0231 static_cast<const G4TouchableHistory*>(aStep->GetPreStepPoint()->GetTouchable());
0232 unsigned int index(0);
0233 int layer(0);
0234 G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
0235
0236 G4ThreeVector localpos = touchable->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
0237 float globalZ = touchable->GetTranslation(0).z();
0238 int iz(globalZ > 0 ? 1 : -1);
0239 int module(-1), cell(-1);
0240 if ((touchable->GetHistoryDepth() == levelT1_) || (touchable->GetHistoryDepth() == levelT2_)) {
0241 layer = touchable->GetReplicaNumber(0);
0242 } else {
0243 layer = touchable->GetReplicaNumber(3);
0244 module = touchable->GetReplicaNumber(2);
0245 cell = touchable->GetReplicaNumber(1);
0246 }
0247 double weight(0);
0248 index = hgcalNumbering_[type]->getUnitID(layer, module, cell, iz, hitPoint, weight);
0249 if (verbosity_ > 1)
0250 edm::LogVerbatim("ValidHGCal") << "HGCal: " << name << " Layer " << layer << " Module " << module << " Cell "
0251 << cell;
0252
0253 double edeposit = aStep->GetTotalEnergyDeposit();
0254 if (verbosity_ > 0)
0255 edm::LogVerbatim("ValidHGCal") << "Layer " << layer << " Index " << std::hex << index << std::dec << " Edep "
0256 << edeposit << " hit at " << hitPoint;
0257 if (detType == 0) {
0258 edepEE_ += edeposit;
0259 if (layer < static_cast<int>(hgcEEedep_.size()))
0260 hgcEEedep_[layer] += edeposit;
0261 } else if (detType == 1) {
0262 edepHEF_ += edeposit;
0263 if (layer < static_cast<int>(hgcHEFedep_.size()))
0264 hgcHEFedep_[layer] += edeposit;
0265 } else {
0266 edepHEB_ += edeposit;
0267 if (layer < static_cast<int>(hgcHEBedep_.size()))
0268 hgcHEBedep_[layer] += edeposit;
0269 }
0270 G4String nextVolume("XXX");
0271 if (aStep->GetTrack()->GetNextVolume() != nullptr)
0272 nextVolume = aStep->GetTrack()->GetNextVolume()->GetName();
0273
0274 if (nextVolume.c_str() != name.c_str()) {
0275 if (std::find(hgchitIndex_.begin(), hgchitIndex_.end(), index) == hgchitIndex_.end()) {
0276 hgchitDets_.emplace_back(dets_[type]);
0277 hgchitIndex_.emplace_back(index);
0278 hgchitX_.emplace_back(hitPoint.x());
0279 hgchitY_.emplace_back(hitPoint.y());
0280 hgchitZ_.emplace_back(hitPoint.z());
0281 }
0282 }
0283 }
0284 }
0285 }
0286 }
0287
0288
0289
0290 void SimG4HGCalValidation::layerAnalysis(PHGCalValidInfo& product) {
0291 if (verbosity_ > 0)
0292 edm::LogVerbatim("ValidHGCal") << "\n ===>>> SimG4HGCalValidation: Energy "
0293 << "deposit\n at EE : " << std::setw(6) << edepEE_ / CLHEP::MeV
0294 << "\n at HEF: " << std::setw(6) << edepHEF_ / CLHEP::MeV
0295 << "\n at HEB: " << std::setw(6) << edepHEB_ / CLHEP::MeV;
0296
0297
0298 product.fillhgcHits(hgchitDets_, hgchitIndex_, hgchitX_, hgchitY_, hgchitZ_);
0299 product.fillhgcLayers(edepEE_, edepHEF_, edepHEB_, hgcEEedep_, hgcHEFedep_, hgcHEBedep_);
0300 }
0301
0302
0303 void SimG4HGCalValidation::clear() {
0304 hgchitDets_.erase(hgchitDets_.begin(), hgchitDets_.end());
0305 hgchitIndex_.erase(hgchitIndex_.begin(), hgchitIndex_.end());
0306 hgchitX_.erase(hgchitX_.begin(), hgchitX_.end());
0307 hgchitY_.erase(hgchitY_.begin(), hgchitY_.end());
0308 hgchitZ_.erase(hgchitZ_.begin(), hgchitZ_.end());
0309 }
0310
0311 DEFINE_SIMWATCHER(SimG4HGCalValidation);