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