Calo/src/HGCSD.cc

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: HGCSD.cc
0003 // Description: Sensitive Detector class for Combined Forward Calorimeter
0004 ///////////////////////////////////////////////////////////////////////////////
0005
0006 #include "DataFormats/Math/interface/FastMath.h"
0007
0008 #include "SimG4CMS/Calo/interface/HGCSD.h"
0009 #include "SimG4Core/Notification/interface/TrackInformation.h"
0010 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0011 #include "FWCore/Utilities/interface/Exception.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0014 #include "Geometry/HGCalTBCommonData/interface/HGCalTBDDDConstants.h"
0015 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
0016 #include "G4LogicalVolumeStore.hh"
0017 #include "G4LogicalVolume.hh"
0018 #include "G4Step.hh"
0019 #include "G4Track.hh"
0020 #include "G4ParticleTable.hh"
0021 #include "G4VProcess.hh"
0022 #include "G4Trap.hh"
0023
0024 #include <fstream>
0025 #include <iomanip>
0026 #include <iostream>
0027 #include <memory>
0028
0029 //#define EDM_ML_DEBUG
0030
0031 HGCSD::HGCSD(const std::string& name,
0032              const HGCalTBDDDConstants* hgc,
0033              const SensitiveDetectorCatalog& clg,
0034              edm::ParameterSet const& p,
0035              const SimTrackManager* manager)
0036     : CaloSD(name,
0037              clg,
0038              p,
0039              manager,
0040              (float)(p.getParameter<edm::ParameterSet>("HGCSD").getParameter<double>("TimeSliceUnit")),
0041              p.getParameter<edm::ParameterSet>("HGCSD").getParameter<bool>("IgnoreTrackID")),
0042       hgcons_(hgc),
0043       slopeMin_(0),
0044       levelT_(99) {
0045 #ifdef plotDebug
0046   tree_ = nullptr;
0047 #endif
0048   numberingScheme_.reset(nullptr);
0049   mouseBite_.reset(nullptr);
0050
0051   edm::ParameterSet m_HGC = p.getParameter<edm::ParameterSet>("HGCSD");
0052   eminHit_ = m_HGC.getParameter<double>("EminHit") * CLHEP::MeV;
0053   storeAllG4Hits_ = m_HGC.getParameter<bool>("StoreAllG4Hits");
0054   rejectMB_ = m_HGC.getParameter<bool>("RejectMouseBite");
0055   waferRot_ = m_HGC.getParameter<bool>("RotatedWafer");
0056   angles_ = m_HGC.getUntrackedParameter<std::vector<double>>("WaferAngles");
0057   double waferSize = m_HGC.getUntrackedParameter<double>("WaferSize") * CLHEP::mm;
0058   double mouseBite = m_HGC.getUntrackedParameter<double>("MouseBite") * CLHEP::mm;
0059   mouseBiteCut_ = waferSize * tan(30.0 * CLHEP::deg) - mouseBite;
0060   dd4hep_ = p.getParameter<bool>("g4GeometryDD4hepSource");
0061
0062   if (storeAllG4Hits_) {
0063     setUseMap(false);
0064     setNumberCheckedHits(0);
0065   }
0066   //this is defined in the hgcsens.xml
0067   G4String myName = name;
0068   myFwdSubdet_ = ForwardSubdetector::ForwardEmpty;
0069   nameX_ = "HGCal";
0070   if (myName.find("HitsEE") != std::string::npos) {
0071     myFwdSubdet_ = ForwardSubdetector::HGCEE;
0072     nameX_ = "HGCalEESensitive";
0073   } else if (myName.find("HitsHEfront") != std::string::npos) {
0074     myFwdSubdet_ = ForwardSubdetector::HGCHEF;
0075     nameX_ = "HGCalHESiliconSensitive";
0076   } else if (myName.find("HitsHEback") != std::string::npos) {
0077     myFwdSubdet_ = ForwardSubdetector::HGCHEB;
0078     nameX_ = "HGCalHEScintillatorSensitive";
0079   }
0080
0081 #ifdef EDM_ML_DEBUG
0082   edm::LogVerbatim("HGCSim") << "**************************************************"
0083                              << "\n"
0084                              << "*                                                *"
0085                              << "\n"
0086                              << "* Constructing a HGCSD  with name " << name << "\n"
0087                              << "*                                                *"
0088                              << "\n"
0089                              << "**************************************************";
0090 #endif
0091   edm::LogVerbatim("HGCSim") << "HGCSD:: Threshold for storing hits: " << eminHit << " for " << nameX_ << " subdet "
0092                              << myFwdSubdet_;
0093   edm::LogVerbatim("HGCSim") << "Flag for storing individual Geant4 Hits " << storeAllG4Hits_;
0094   edm::LogVerbatim("HGCSim") << "Reject MosueBite Flag: " << rejectMB_ << " Size of wafer " << waferSize
0095                              << " Mouse Bite " << mouseBite << ":" << mouseBiteCut_ << " along " << angles_.size()
0096                              << " axes";
0097 }
0098
0099 double HGCSD::getEnergyDeposit(const G4Step* aStep) {
0100   double r = aStep->GetPreStepPoint()->GetPosition().perp();
0101   double z = std::abs(aStep->GetPreStepPoint()->GetPosition().z());
0102
0103 #ifdef EDM_ML_DEBUG
0104   G4int parCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
0105   G4LogicalVolume* lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
0106   edm::LogVerbatim("HGCSim") << "HGCSD: Hit from standard path from " << lv->GetName() << " for Track "
0107                              << aStep->GetTrack()->GetTrackID() << " ("
0108                              << aStep->GetTrack()->GetDefinition()->GetParticleName() << ":" << parCode << ") R = " << r
0109                              << " Z = " << z << " slope = " << r / z << ":" << slopeMin_;
0110 #endif
0111
0112   // Apply fiductial volume
0113   if (r < z * slopeMin_) {
0114     return 0.0;
0115   }
0116
0117   double wt1 = getResponseWt(aStep->GetTrack());
0118   double wt2 = aStep->GetTrack()->GetWeight();
0119   double destep = wt1 * aStep->GetTotalEnergyDeposit();
0120   if (wt2 > 0)
0121     destep *= wt2;
0122
0123 #ifdef plotDebug
0124   const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0125   G4double tmptrackE = aStep->GetTrack()->GetKineticEnergy();
0126   G4int parCodex = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
0127   G4double angle = (aStep->GetTrack()->GetMomentumDirection().theta()) / CLHEP::deg;
0128   G4int layer = ((touch->GetHistoryDepth() == levelT_) ? touch->GetReplicaNumber(0) : touch->GetReplicaNumber(2));
0129   G4int ilayer = (layer - 1) / 3;
0130   if (aStep->GetTotalEnergyDeposit() > 0) {
0131     t_Layer_.emplace_back(ilayer);
0132     t_Parcode_.emplace_back(parCodex);
0133     t_dEStep1_.emplace_back(aStep->GetTotalEnergyDeposit());
0134     t_dEStep2_.emplace_back(destep);
0135     t_TrackE_.emplace_back(tmptrackE);
0136     t_Angle_.emplace_back(angle);
0137   }
0138 #endif
0139
0140   return destep;
0141 }
0142
0143 uint32_t HGCSD::setDetUnitId(const G4Step* aStep) {
0144   const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0145   const G4VTouchable* touch = preStepPoint->GetTouchable();
0146
0147   //determine the exact position in global coordinates in the mass geometry
0148   G4ThreeVector hitPoint = preStepPoint->GetPosition();
0149   float globalZ = touch->GetTranslation(0).z();
0150   int iz(globalZ > 0 ? 1 : -1);
0151
0152   //convert to local coordinates (=local to the current volume):
0153   G4ThreeVector localpos = touch->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
0154
0155   //get the det unit id with
0156   ForwardSubdetector subdet = myFwdSubdet_;
0157
0158   int layer(-1), moduleLev(-1), module(-1), cell(-1);
0159   if (touch->GetHistoryDepth() == levelT_) {
0160     layer = touch->GetReplicaNumber(0);
0161 #ifdef EDM_ML_DEBUG
0162     edm::LogVerbatim("HGCSim") << "Depths: " << touch->GetHistoryDepth() << " name " << touch->GetVolume(0)->GetName()
0163                                << " layer:module:cell " << layer << ":" << moduleLev << ":" << module << ":" << cell;
0164 #endif
0165   } else {
0166     layer = touch->GetReplicaNumber(2);
0167     module = touch->GetReplicaNumber(1);
0168     cell = touch->GetReplicaNumber(0);
0169     moduleLev = 1;
0170   }
0171 #ifdef EDM_ML_DEBUG
0172   const G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
0173   edm::LogVerbatim("HGCSim") << "Depths: " << touch->GetHistoryDepth() << " name " << touch->GetVolume(0)->GetName()
0174                              << ":" << touch->GetReplicaNumber(0) << "   " << touch->GetVolume(1)->GetName() << ":"
0175                              << touch->GetReplicaNumber(1) << "   " << touch->GetVolume(2)->GetName() << ":"
0176                              << touch->GetReplicaNumber(2) << "    layer:module:cell " << layer << ":" << moduleLev
0177                              << ":" << module << ":" << cell << " Material " << mat->GetName() << ":"
0178                              << mat->GetRadlen();
0179   for (int k = 0; k < touch->GetHistoryDepth(); ++k)
0180     edm::LogVerbatim("HGCSim") << "Level [" << k << "] " << touch->GetVolume(k)->GetName() << ":"
0181                                << touch->GetReplicaNumber(k);
0182 #endif
0183   // The following statement should be examined later before elimination
0184   // VI: this is likely a check if media is vacuum - not needed
0185   if (aStep->GetPreStepPoint()->GetMaterial()->GetRadlen() > 100000.)
0186     return 0;
0187
0188   uint32_t id = setDetUnitId(subdet, layer, module, cell, iz, localpos);
0189   if (rejectMB_ && id != 0) {
0190     int det, z, lay, wafer, type, ic;
0191     HGCalTestNumbering::unpackHexagonIndex(id, det, z, lay, wafer, type, ic);
0192 #ifdef EDM_ML_DEBUG
0193     edm::LogVerbatim("HGCSim") << "ID " << std::hex << id << std::dec << " Input " << subdet << ":" << layer << ":"
0194                                << module << ":" << cell << ":" << iz << localpos.x() << ":" << localpos.y()
0195                                << " Decode " << det << ":" << z << ":" << lay << ":" << wafer << ":" << type << ":"
0196                                << ic;
0197 #endif
0198     G4ThreeVector local =
0199         ((moduleLev >= 0) ? (touch->GetHistory()->GetTransform(moduleLev).TransformPoint(hitPoint)) : G4ThreeVector());
0200     if (mouseBite_->exclude(local, z, layer, wafer, 0))
0201       id = 0;
0202   }
0203   return id;
0204 }
0205
0206 void HGCSD::update(const BeginOfJob* job) {
0207   if (hgcons_ != nullptr) {
0208     geom_mode_ = hgcons_->geomMode();
0209     slopeMin_ = hgcons_->minSlope();
0210     levelT_ = hgcons_->levelTop();
0211     if (dd4hep_)
0212       ++levelT_;
0213     numberingScheme_ = std::make_unique<HGCNumberingScheme>(*hgcons_, nameX_);
0214     if (rejectMB_)
0215       mouseBite_ = std::make_unique<HGCMouseBite>(*hgcons_, angles_, mouseBiteCut_, waferRot_);
0216   } else {
0217     edm::LogError("HGCSim") << "HGCSD : Cannot find HGCalTBDDDConstants for " << nameX_;
0218     throw cms::Exception("Unknown", "HGCSD") << "Cannot find HGCalTBDDDConstants for " << nameX_ << "\n";
0219   }
0220 #ifdef EDM_ML_DEBUG
0221   edm::LogVerbatim("HGCSim") << "HGCSD::Initialized with mode " << geom_mode_ << " Slope cut " << slopeMin_
0222                              << " top Level " << levelT_ << " dd4hep " << dd4hep_;
0223 #endif
0224 }
0225
0226 void HGCSD::initRun() {
0227 #ifdef plotDebug
0228   edm::Service<TFileService> tfile;
0229   if (tfile.isAvailable()) {
0230     tree_ = tfile->make<TTree>("TreeHGCSD", "TreeHGCSD");
0231     tree_->Branch("EventID", &t_EventID_);
0232     tree_->Branch("Layer", &t_Layer_);
0233     tree_->Branch("ParticleCode", &t_Parcode_);
0234     tree_->Branch("dEStepOriginal", &t_dEStep1_);
0235     tree_->Branch("dEStepWeighted", &t_dEStep2_);
0236     tree_->Branch("TrackEnergy", &t_TrackE_);
0237     tree_->Branch("ThetaAngle", &t_Angle_);
0238   }
0239 #endif
0240 }
0241
0242 void HGCSD::initEvent(const BeginOfEvent* g4Event) {
0243   const G4Event* evt = (*g4Event)();
0244   t_EventID_ = evt->GetEventID();
0245 #ifdef plotDebug
0246   t_Layer_.clear();
0247   t_Parcode_.clear();
0248   t_dEStep1_.clear();
0249   t_dEStep2_.clear();
0250   t_TrackE_.clear();
0251   t_Angle_.clear();
0252 #endif
0253 }
0254
0255 void HGCSD::endEvent() {
0256 #ifdef plotDebug
0257   if (tree_)
0258     tree_->Fill();
0259 #endif
0260 }
0261
0262 bool HGCSD::filterHit(CaloG4Hit* aHit, double time) {
0263   return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > eminHit_));
0264 }
0265
0266 uint32_t HGCSD::setDetUnitId(ForwardSubdetector& subdet, int layer, int module, int cell, int iz, G4ThreeVector& pos) {
0267   uint32_t id = numberingScheme_ ? numberingScheme_->getUnitID(subdet, layer, module, cell, iz, pos) : 0;
0268   return id;
0269 }