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/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 // to retreive hits
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&);  // stop default
0062   const SimG4HGCalValidation& operator=(const SimG4HGCalValidation&);
0063
0064   void init();
0065
0066   // observer classes
0067   void update(const BeginOfEvent* evt) override;
0068   void update(const G4Step* step) override;
0069
0070   // analysis related class
0071   void layerAnalysis(PHGCalValidInfo&);
0072   void clear();
0073
0074 private:
0075   //HGCal numbering scheme
0076   std::vector<edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> > ddconsToken_;
0077   std::vector<HGCNumberingScheme*> hgcNumbering_;
0078   std::vector<HGCalNumberingScheme*> hgcalNumbering_;
0079
0080   // to read from parameter set
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   // parameters from geometry
0087   int levelT1_, levelT2_;
0088
0089   // some private members for ananlysis
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 //=================================================================== per EVENT
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   //HGCal variables
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   //Cache reset
0230   clear();
0231 }
0232
0233 //=================================================================== each STEP
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     // Only for Sensitive detector
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       // Right type of SD
0257       if (type >= 0) {
0258         //Get the 32-bit index of the hit cell
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         // HGCal
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()) {  //save hit when it exits cell
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       }  // it is right type of SD
0325     }    // it is in a SD
0326   }      //if aStep!=NULL
0327 }  //end update aStep
0328
0329 //================================================================ End of EVENT
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   //Fill HGC Variables
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);