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