HGCalValidation/plugins/SimG4HGCalValidation.cc

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: SimG4HGCalValidation.cc
0003 // Description: Main analysis class for HGCal Validation of G4 Hits
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 // to retreive hits
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&);  // stop default
0055   const SimG4HGCalValidation& operator=(const SimG4HGCalValidation&);
0056
0057   void init();
0058
0059   // observer classes
0060   void update(const BeginOfEvent* evt) override;
0061   void update(const G4Step* step) override;
0062
0063   // analysis related class
0064   void layerAnalysis(PHGCalValidInfo&);
0065   void clear();
0066
0067 private:
0068   //HGCal numbering scheme
0069   std::vector<edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord>> ddconsToken_;
0070   std::vector<std::unique_ptr<HGCalNumberingScheme>> hgcalNumbering_;
0071
0072   // to read from parameter set
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   // parameters from geometry
0082   int levelT1_, levelT2_;
0083
0084   // some private members for ananlysis
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 //=================================================================== per EVENT
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   //HGCal variables
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   //Cache reset
0201   clear();
0202 }
0203
0204 //=================================================================== each STEP
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     // Only for Sensitive detector
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       // Right type of SD
0228       if (type >= 0) {
0229         //Get the 32-bit index of the hit cell
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         // HGCal
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()) {  //save hit when it exits cell
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       }  // it is right type of SD
0284     }  // it is in a SD
0285   }  //if aStep!=NULL
0286 }  //end update aStep
0287
0288 //================================================================ End of EVENT
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   //Fill HGC Variables
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);