Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-06-29 02:26:08 

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: HGCalSD.cc
0003 // Description: Sensitive Detector class for High Granularity Calorimeter
0004 ///////////////////////////////////////////////////////////////////////////////
0005 
0006 #include "DataFormats/Math/interface/FastMath.h"
0007 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0008 #include "SimG4CMS/Calo/interface/HGCalSD.h"
0009 #include "SimG4Core/Notification/interface/TrackInformation.h"
0010 #include "FWCore/Utilities/interface/Exception.h"
0011 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0012 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
0013 #include "G4LogicalVolumeStore.hh"
0014 #include "G4LogicalVolume.hh"
0015 #include "G4Step.hh"
0016 #include "G4Track.hh"
0017 #include "G4ParticleTable.hh"
0018 #include "G4VProcess.hh"
0019 #include "G4Trap.hh"
0020 
0021 #include <fstream>
0022 #include <iomanip>
0023 #include <iostream>
0024 #include <memory>
0025 
0026 //#define EDM_ML_DEBUG
0027 
0028 HGCalSD::HGCalSD(const std::string& name,
0029                  const HGCalDDDConstants* hgc,
0030                  const SensitiveDetectorCatalog& clg,
0031                  edm::ParameterSet const& p,
0032                  const SimTrackManager* manager)
0033     : CaloSD(name,
0034              clg,
0035              p,
0036              manager,
0037              static_cast<float>(p.getParameter<edm::ParameterSet>("HGCSD").getParameter<double>("TimeSliceUnit")),
0038              p.getParameter<edm::ParameterSet>("HGCSD").getParameter<bool>("IgnoreTrackID")),
0039       hgcons_(hgc),
0040       slopeMin_(0),
0041       levelT1_(99),
0042       levelT2_(99),
0043       useSimWt_(0),
0044       tan30deg_(std::tan(30.0 * CLHEP::deg)) {
0045   numberingScheme_.reset(nullptr);
0046   mouseBite_.reset(nullptr);
0047 
0048   edm::ParameterSet m_HGC = p.getParameter<edm::ParameterSet>("HGCSD");
0049   eminHit_ = m_HGC.getParameter<double>("EminHit") * CLHEP::MeV;
0050   fiducialCut_ = m_HGC.getParameter<bool>("FiducialCut");
0051   distanceFromEdge_ = m_HGC.getParameter<double>("DistanceFromEdge");
0052   storeAllG4Hits_ = m_HGC.getParameter<bool>("StoreAllG4Hits");
0053   rejectMB_ = m_HGC.getParameter<bool>("RejectMouseBite");
0054   waferRot_ = m_HGC.getParameter<bool>("RotatedWafer");
0055   cornerMinMask_ = m_HGC.getParameter<int>("CornerMinMask");
0056   angles_ = m_HGC.getUntrackedParameter<std::vector<double>>("WaferAngles");
0057 
0058   if (storeAllG4Hits_) {
0059     setUseMap(false);
0060     setNumberCheckedHits(0);
0061   }
0062 
0063   //this is defined in the hgcsens.xml
0064   G4String myName = name;
0065   mydet_ = DetId::Forward;
0066   nameX_ = "HGCal";
0067   if (myName.find("HitsEE") != std::string::npos) {
0068     mydet_ = DetId::HGCalEE;
0069     nameX_ = "HGCalEESensitive";
0070   } else if (myName.find("HitsHEfront") != std::string::npos) {
0071     mydet_ = DetId::HGCalHSi;
0072     nameX_ = "HGCalHESiliconSensitive";
0073   }
0074 
0075 #ifdef EDM_ML_DEBUG
0076   edm::LogVerbatim("HGCSim") << "**************************************************"
0077                              << "\n"
0078                              << "*                                                *"
0079                              << "\n"
0080                              << "* Constructing a HGCalSD  with name " << name << "\n"
0081                              << "*                                                *"
0082                              << "\n"
0083                              << "**************************************************";
0084 #endif
0085   edm::LogVerbatim("HGCSim") << "HGCalSD:: Threshold for storing hits: " << eminHit_ << " for " << nameX_
0086                              << " detector " << mydet_;
0087   edm::LogVerbatim("HGCSim") << "Flag for storing individual Geant4 Hits " << storeAllG4Hits_;
0088   edm::LogVerbatim("HGCSim") << "Fiducial volume cut with cut from eta/phi "
0089                              << "boundary " << fiducialCut_ << " at " << distanceFromEdge_;
0090   edm::LogVerbatim("HGCSim") << "Reject MosueBite Flag: " << rejectMB_ << " cuts along " << angles_.size()
0091                              << " axes: " << angles_[0] << ", " << angles_[1];
0092 }
0093 
0094 double HGCalSD::getEnergyDeposit(const G4Step* aStep) {
0095   double r = aStep->GetPreStepPoint()->GetPosition().perp();
0096   double z = std::abs(aStep->GetPreStepPoint()->GetPosition().z());
0097 #ifdef EDM_ML_DEBUG
0098   G4int parCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
0099   G4String parName = aStep->GetTrack()->GetDefinition()->GetParticleName();
0100   G4LogicalVolume* lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
0101   edm::LogVerbatim("HGCSim") << "HGCalSD: Hit from standard path from " << lv->GetName() << " for Track "
0102                              << aStep->GetTrack()->GetTrackID() << " (" << parCode << ":" << parName << ") R = " << r
0103                              << " Z = " << z << " slope = " << r / z << ":" << slopeMin_;
0104 #endif
0105   // Apply fiducial cut
0106   if (r < z * slopeMin_) {
0107 #ifdef EDM_ML_DEBUG
0108     edm::LogVerbatim("HGCSim") << "HGCalSD: Fiducial Volume cut";
0109 #endif
0110     return 0.0;
0111   }
0112 
0113   double wt1 = getResponseWt(aStep->GetTrack());
0114   double wt2 = aStep->GetTrack()->GetWeight();
0115   double destep = weight_ * wt1 * (aStep->GetTotalEnergyDeposit());
0116   if (wt2 > 0)
0117     destep *= wt2;
0118 #ifdef EDM_ML_DEBUG
0119   edm::LogVerbatim("HGCSim") << "HGCalSD: weights= " << weight_ << ":" << wt1 << ":" << wt2 << " Total weight "
0120                              << weight_ * wt1 * wt2 << " deStep: " << aStep->GetTotalEnergyDeposit() << ":" << destep;
0121 #endif
0122   return destep;
0123 }
0124 
0125 uint32_t HGCalSD::setDetUnitId(const G4Step* aStep) {
0126   const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0127   const G4VTouchable* touch = preStepPoint->GetTouchable();
0128 
0129 #ifdef EDM_ML_DEBUG
0130   edm::LogVerbatim("HGCSim") << "DepthsTop: " << touch->GetHistoryDepth() << ":" << levelT1_ << ":" << levelT2_;
0131   printDetectorLevels(touch);
0132 #endif
0133   //determine the exact position in global coordinates in the mass geometry
0134   G4ThreeVector hitPoint = preStepPoint->GetPosition();
0135   float globalZ = touch->GetTranslation(0).z();
0136   int iz(globalZ > 0 ? 1 : -1);
0137 
0138   int layer(0), module(-1), cell(-1);
0139   if ((geom_mode_ == HGCalGeometryMode::Hexagon8Module) || (geom_mode_ == HGCalGeometryMode::Hexagon8Cassette)) {
0140     if (useSimWt_ > 0) {
0141       layer = touch->GetReplicaNumber(2);
0142       module = touch->GetReplicaNumber(1);
0143     } else if (touch->GetHistoryDepth() > levelT2_) {
0144       layer = touch->GetReplicaNumber(4);
0145       module = touch->GetReplicaNumber(3);
0146       cell = touch->GetReplicaNumber(1);
0147     } else {
0148       layer = touch->GetReplicaNumber(3);
0149       module = touch->GetReplicaNumber(2);
0150     }
0151 #ifdef EDM_ML_DEBUG
0152     edm::LogVerbatim("HGCSim") << "DepthsTop: " << touch->GetHistoryDepth() << ":" << levelT1_ << ":" << levelT2_ << ":"
0153                                << useSimWt_ << " name " << touch->GetVolume(0)->GetName() << " layer:module:cell "
0154                                << layer << ":" << module << ":" << cell;
0155     printDetectorLevels(touch);
0156 #endif
0157   } else if ((touch->GetHistoryDepth() == levelT1_) || (touch->GetHistoryDepth() == levelT2_)) {
0158     layer = touch->GetReplicaNumber(0);
0159 #ifdef EDM_ML_DEBUG
0160     edm::LogVerbatim("HGCSim") << "DepthsTop: " << touch->GetHistoryDepth() << ":" << levelT1_ << ":" << levelT2_
0161                                << " name " << touch->GetVolume(0)->GetName() << " layer:module:cell " << layer << ":"
0162                                << module << ":" << cell;
0163 #endif
0164   } else {
0165     layer = touch->GetReplicaNumber(3);
0166     module = touch->GetReplicaNumber(2);
0167     cell = touch->GetReplicaNumber(1);
0168 #ifdef EDM_ML_DEBUG
0169     edm::LogVerbatim("HGCSim") << "DepthsInside: " << touch->GetHistoryDepth() << " name "
0170                                << touch->GetVolume(0)->GetName() << " layer:module:cell " << layer << ":" << module
0171                                << ":" << cell;
0172 #endif
0173   }
0174 #ifdef EDM_ML_DEBUG
0175   G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
0176   edm::LogVerbatim("HGCSim") << "Depths: " << touch->GetHistoryDepth() << " name " << touch->GetVolume(0)->GetName()
0177                              << ":" << touch->GetReplicaNumber(0) << "   " << touch->GetVolume(1)->GetName() << ":"
0178                              << touch->GetReplicaNumber(1) << "   " << touch->GetVolume(2)->GetName() << ":"
0179                              << touch->GetReplicaNumber(2) << "   " << touch->GetVolume(3)->GetName() << ":"
0180                              << touch->GetReplicaNumber(3) << "   " << touch->GetVolume(4)->GetName() << ":"
0181                              << touch->GetReplicaNumber(4) << "   "
0182                              << " layer:module:cell " << layer << ":" << module << ":" << cell << " Material "
0183                              << mat->GetName() << ":" << mat->GetRadlen();
0184 #endif
0185   // The following statement should be examined later before elimination
0186   if (aStep->GetPreStepPoint()->GetMaterial()->GetRadlen() > 100000.)
0187     return 0;
0188 
0189   uint32_t id = setDetUnitId(layer, module, cell, iz, hitPoint);
0190   if (rejectMB_ && id != 0) {
0191     auto uv = HGCSiliconDetId(id).waferUV();
0192 #ifdef EDM_ML_DEBUG
0193     edm::LogVerbatim("HGCSim") << "ID " << std::hex << id << std::dec << " " << HGCSiliconDetId(id);
0194 #endif
0195     if (mouseBite_->exclude(hitPoint, iz, uv.first, uv.second)) {
0196       id = 0;
0197 #ifdef EDM_ML_DEBUG
0198       edm::LogVerbatim("HGCSim") << "Rejected by mousebite cutoff *****";
0199 #endif
0200     }
0201   }
0202 #ifdef EDM_ML_DEBUG
0203   if (id != 0)
0204     edm::LogVerbatim("HGCSim") << HGCSiliconDetId(id);
0205 #endif
0206   return id;
0207 }
0208 
0209 void HGCalSD::update(const BeginOfJob* job) {
0210   if (hgcons_ != nullptr) {
0211     geom_mode_ = hgcons_->geomMode();
0212     slopeMin_ = hgcons_->minSlope();
0213     levelT1_ = hgcons_->levelTop(0);
0214     levelT2_ = hgcons_->levelTop(1);
0215     useSimWt_ = hgcons_->getParameter()->useSimWt_;
0216     double waferSize = hgcons_->waferSize(false);
0217     double mouseBite = hgcons_->mouseBite(false);
0218     mouseBiteCut_ = waferSize * tan30deg_ - mouseBite;
0219 #ifdef EDM_ML_DEBUG
0220     edm::LogVerbatim("HGCSim") << "HGCalSD::Initialized with mode " << geom_mode_ << " Slope cut " << slopeMin_
0221                                << " top Level " << levelT1_ << ":" << levelT2_ << " useSimWt " << useSimWt_ << " wafer "
0222                                << waferSize << ":" << mouseBite;
0223 #endif
0224 
0225     numberingScheme_ = std::make_unique<HGCalNumberingScheme>(*hgcons_, mydet_, nameX_);
0226     if (rejectMB_)
0227       mouseBite_ = std::make_unique<HGCMouseBite>(*hgcons_, angles_, mouseBiteCut_, waferRot_);
0228   } else {
0229     throw cms::Exception("Unknown", "HGCalSD") << "Cannot find HGCalDDDConstants for " << nameX_ << "\n";
0230   }
0231 }
0232 
0233 void HGCalSD::initRun() {}
0234 
0235 bool HGCalSD::filterHit(CaloG4Hit* aHit, double time) {
0236   return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > eminHit_));
0237 }
0238 
0239 uint32_t HGCalSD::setDetUnitId(int layer, int module, int cell, int iz, G4ThreeVector& pos) {
0240   uint32_t id = numberingScheme_ ? numberingScheme_->getUnitID(layer, module, cell, iz, pos, weight_) : 0;
0241   if (cornerMinMask_ > 2) {
0242     if (hgcons_->maskCell(DetId(id), cornerMinMask_)) {
0243       id = 0;
0244       ignoreRejection();
0245     }
0246   }
0247   if (hgcons_->waferHexagon8File() || (id == 0))
0248     ignoreRejection();
0249   return id;
0250 }
0251 
0252 bool HGCalSD::isItinFidVolume(const G4ThreeVector& pos) {
0253   if (fiducialCut_) {
0254     return (hgcons_->distFromEdgeHex(pos.x(), pos.y(), pos.z()) > distanceFromEdge_);
0255   } else {
0256     return true;
0257   }
0258 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!