Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:29:56 

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: HGCScintSD.cc
0003 // Description: Sensitive Detector class for the Scintillator part of
0004 //              High Granularity Calorimeter
0005 ///////////////////////////////////////////////////////////////////////////////
0006 
0007 #include "DataFormats/Math/interface/FastMath.h"
0008 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0009 #include "SimG4CMS/Calo/interface/HGCScintSD.h"
0010 #include "SimG4CMS/Calo/interface/CaloSimUtils.h"
0011 #include "SimG4Core/Notification/interface/TrackInformation.h"
0012 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0013 #include "FWCore/ParameterSet/interface/FileInPath.h"
0014 #include "FWCore/Utilities/interface/Exception.h"
0015 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0016 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
0017 #include "G4LogicalVolumeStore.hh"
0018 #include "G4LogicalVolume.hh"
0019 #include "G4Step.hh"
0020 #include "G4Track.hh"
0021 #include "G4ParticleTable.hh"
0022 #include "G4VProcess.hh"
0023 #include "G4Trap.hh"
0024 
0025 #include <fstream>
0026 #include <iomanip>
0027 #include <iostream>
0028 #include <memory>
0029 
0030 //#define EDM_ML_DEBUG
0031 
0032 HGCScintSD::HGCScintSD(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       levelT1_(99),
0046       levelT2_(99),
0047       firstLayer_(0) {
0048   numberingScheme_.reset(nullptr);
0049 
0050   edm::ParameterSet m_HGC = p.getParameter<edm::ParameterSet>("HGCScintSD");
0051   eminHit_ = m_HGC.getParameter<double>("EminHit") * CLHEP::MeV;
0052   fiducialCut_ = m_HGC.getParameter<bool>("FiducialCut");
0053   distanceFromEdge_ = m_HGC.getParameter<double>("DistanceFromEdge");
0054   useBirk_ = m_HGC.getParameter<bool>("UseBirkLaw");
0055   birk1_ = m_HGC.getParameter<double>("BirkC1") * (CLHEP::g / (CLHEP::MeV * CLHEP::cm2));
0056   birk2_ = m_HGC.getParameter<double>("BirkC2");
0057   birk3_ = m_HGC.getParameter<double>("BirkC3");
0058   storeAllG4Hits_ = m_HGC.getParameter<bool>("StoreAllG4Hits");
0059   checkID_ = m_HGC.getUntrackedParameter<bool>("CheckID");
0060   fileName_ = m_HGC.getUntrackedParameter<std::string>("TileFileName");
0061   verbose_ = m_HGC.getUntrackedParameter<int>("Verbosity");
0062 
0063   if (storeAllG4Hits_) {
0064     setUseMap(false);
0065     setNumberCheckedHits(0);
0066   }
0067 
0068   //this is defined in the hgcsens.xml
0069   G4String myName = name;
0070   mydet_ = DetId::Forward;
0071   nameX_ = "HGCal";
0072   if (myName.find("HitsHEback") != std::string::npos) {
0073     mydet_ = DetId::HGCalHSc;
0074     nameX_ = "HGCalHEScintillatorSensitive";
0075   }
0076 
0077 #ifdef EDM_ML_DEBUG
0078   edm::LogVerbatim("HGCSim") << "**************************************************"
0079                              << "\n"
0080                              << "*                                                *"
0081                              << "\n"
0082                              << "* Constructing a HGCScintSD  with name " << name << "\n"
0083                              << "*                                                *"
0084                              << "\n"
0085                              << "**************************************************";
0086 #endif
0087   edm::LogVerbatim("HGCSim") << "HGCScintSD:: Threshold for storing hits: " << eminHit_ << " for " << nameX_
0088                              << " detector " << mydet_ << " File " << fileName_;
0089   edm::LogVerbatim("HGCSim") << "Flag for storing individual Geant4 Hits " << storeAllG4Hits_;
0090   edm::LogVerbatim("HGCSim") << "Fiducial volume cut with cut from eta/phi "
0091                              << "boundary " << fiducialCut_ << " at " << distanceFromEdge_;
0092   edm::LogVerbatim("HGCSim") << "Use of Birks law is set to      " << useBirk_
0093                              << "  with three constants kB = " << birk1_ << ", C1 = " << birk2_ << ", C2 = " << birk3_;
0094 
0095   if (!fileName_.empty()) {
0096     edm::FileInPath filetmp("SimG4CMS/Calo/data/" + fileName_);
0097     std::string fileName = filetmp.fullPath();
0098     std::ifstream fInput(fileName.c_str());
0099     if (!fInput.good()) {
0100       edm::LogVerbatim("HGCSim") << "Cannot open file " << fileName;
0101     } else {
0102       char buffer[80];
0103       while (fInput.getline(buffer, 80)) {
0104         std::vector<std::string> items = CaloSimUtils::splitString(std::string(buffer));
0105         if (items.size() > 2) {
0106           int layer = std::atoi(items[0].c_str());
0107           int ring = std::atoi(items[1].c_str());
0108           int phi = std::atoi(items[2].c_str());
0109           tiles_.emplace_back(HGCalTileIndex::tileIndex(layer, ring, phi));
0110         }
0111       }
0112       edm::LogVerbatim("HGCSim") << "Reads in " << tiles_.size() << " tile information from " << fileName_;
0113       fInput.close();
0114     }
0115   }
0116 }
0117 
0118 double HGCScintSD::getEnergyDeposit(const G4Step* aStep) {
0119   double r = aStep->GetPreStepPoint()->GetPosition().perp();
0120   double z = std::abs(aStep->GetPreStepPoint()->GetPosition().z());
0121 #ifdef EDM_ML_DEBUG
0122   G4int parCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
0123   G4String parName = aStep->GetTrack()->GetDefinition()->GetParticleName();
0124   G4LogicalVolume* lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
0125   edm::LogVerbatim("HGCSim") << "HGCScintSD: Hit from standard path from " << lv->GetName() << " for Track "
0126                              << aStep->GetTrack()->GetTrackID() << " (" << parCode << ":" << parName << ") R = " << r
0127                              << " Z = " << z << " slope = " << r / z << ":" << slopeMin_;
0128 #endif
0129   // Apply fiducial cut
0130   if (r < z * slopeMin_) {
0131 #ifdef EDM_ML_DEBUG
0132     edm::LogVerbatim("HGCSim") << "HGCScintSD: Fiducial Volume cut";
0133 #endif
0134     return 0.0;
0135   }
0136 
0137   double wt1 = getResponseWt(aStep->GetTrack());
0138   double wt2 = aStep->GetTrack()->GetWeight();
0139   double wt3 = (useBirk_ ? getAttenuation(aStep, birk1_, birk2_, birk3_) : 1.0);
0140   double destep = weight_ * wt1 * wt3 * (aStep->GetTotalEnergyDeposit());
0141   if (wt2 > 0)
0142     destep *= wt2;
0143 #ifdef EDM_ML_DEBUG
0144   edm::LogVerbatim("HGCalSim") << "HGCScintSD: weights= " << weight_ << ":" << wt1 << ":" << wt2 << ":" << wt3
0145                                << " Total weight " << weight_ * wt1 * wt2 * wt3
0146                                << " deStep: " << aStep->GetTotalEnergyDeposit() << ":" << destep;
0147 #endif
0148   return destep;
0149 }
0150 
0151 uint32_t HGCScintSD::setDetUnitId(const G4Step* aStep) {
0152   const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0153   const G4VTouchable* touch = preStepPoint->GetTouchable();
0154 
0155 #ifdef EDM_ML_DEBUG
0156   edm::LogVerbatim("HGCSim") << "DepthsTop: " << touch->GetHistoryDepth() << ":" << levelT1_ << ":" << levelT2_;
0157   printDetectorLevels(touch);
0158 #endif
0159   //determine the exact position in global coordinates in the mass geometry
0160   G4ThreeVector hitPoint = preStepPoint->GetPosition();
0161   float globalZ = touch->GetTranslation(0).z();
0162   int iz(globalZ > 0 ? 1 : -1);
0163 
0164   int layer(0), module(-1), cell(-1);
0165   if ((geom_mode_ == HGCalGeometryMode::TrapezoidModule) || (geom_mode_ == HGCalGeometryMode::TrapezoidCassette)) {
0166     layer = touch->GetReplicaNumber(1);
0167   } else if ((touch->GetHistoryDepth() == levelT1_) || (touch->GetHistoryDepth() == levelT2_)) {
0168     layer = touch->GetReplicaNumber(0);
0169 #ifdef EDM_ML_DEBUG
0170     edm::LogVerbatim("HGCSim") << "DepthsTop: " << touch->GetHistoryDepth() << ":" << levelT1_ << ":" << levelT2_
0171                                << " name " << touch->GetVolume(0)->GetName() << " layer:module:cell " << layer << ":"
0172                                << module << ":" << cell;
0173 #endif
0174   } else {
0175     layer = touch->GetReplicaNumber(3);
0176     module = touch->GetReplicaNumber(2);
0177     cell = touch->GetReplicaNumber(1);
0178 #ifdef EDM_ML_DEBUG
0179     edm::LogVerbatim("HGCSim") << "DepthsInside: " << touch->GetHistoryDepth() << " name "
0180                                << touch->GetVolume(0)->GetName() << " layer:module:cell " << layer << ":" << module
0181                                << ":" << cell;
0182 #endif
0183   }
0184 #ifdef EDM_ML_DEBUG
0185   G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
0186   edm::LogVerbatim("HGCSim") << "Depths: " << touch->GetHistoryDepth() << " name " << touch->GetVolume(0)->GetName()
0187                              << ":" << touch->GetReplicaNumber(0) << "   " << touch->GetVolume(1)->GetName() << ":"
0188                              << touch->GetReplicaNumber(1) << "   " << touch->GetVolume(2)->GetName() << ":"
0189                              << touch->GetReplicaNumber(2) << "   " << touch->GetVolume(3)->GetName() << ":"
0190                              << touch->GetReplicaNumber(3) << "   " << touch->GetVolume(4)->GetName() << ":"
0191                              << touch->GetReplicaNumber(4) << "   "
0192                              << " layer:module:cell " << layer << ":" << module << ":" << cell << " Material "
0193                              << mat->GetName() << ":" << mat->GetRadlen();
0194 #endif
0195   // The following statement should be examined later before elimination
0196   if (aStep->GetPreStepPoint()->GetMaterial()->GetRadlen() > 100000.)
0197     return 0;
0198 
0199   uint32_t id = setDetUnitId(layer, module, cell, iz, hitPoint);
0200   bool debug(false);
0201   if (!tiles_.empty()) {
0202     HGCScintillatorDetId hid(id);
0203     int indx = HGCalTileIndex::tileIndex(firstLayer_ + hid.layer(), hid.ring(), hid.iphi());
0204     if (std::find(tiles_.begin(), tiles_.end(), indx) != tiles_.end())
0205       debug = true;
0206   }
0207   if (debug)
0208     edm::LogVerbatim("HGCSim") << "Layer:module:cell:iz " << layer << ":" << module << ":" << cell << ":" << iz
0209                                << "  Point (" << hitPoint.x() << ", " << hitPoint.y() << ", " << hitPoint.z() << ") "
0210                                << HGCScintillatorDetId(id);
0211 
0212   if (!isItinFidVolume(hitPoint)) {
0213 #ifdef EDM_ML_DEBUG
0214     edm::LogVerbatim("HGCSim") << "ID " << std::hex << id << std::dec << " " << HGCScintillatorDetId(id)
0215                                << " is rejected by fiducilal volume cut";
0216 #endif
0217     id = 0;
0218   }
0219   if ((id != 0) && checkID_) {
0220     HGCScintillatorDetId hid1(id);
0221     std::string_view pid =
0222         ((hgcons_->cassetteShiftScintillator(hid1.zside(), hid1.layer(), hid1.iphi())) ? "HGCSim" : "HGCalSim");
0223     bool debug = (verbose_ > 0) ? true : false;
0224     auto xy = hgcons_->locateCell(HGCScintillatorDetId(id), debug);
0225     double dx = xy.first - (hitPoint.x() / CLHEP::cm);
0226     double dy = xy.second - (hitPoint.y() / CLHEP::cm);
0227     double diff = std::sqrt(dx * dx + dy * dy);
0228     constexpr double tol = 10.0;
0229     bool valid = hgcons_->isValidTrap(hid1.zside(), hid1.layer(), hid1.ring(), hid1.iphi());
0230     if ((diff > tol) || (!valid))
0231       pid = "HGCalError";
0232     edm::LogVerbatim(pid) << "CheckID " << HGCScintillatorDetId(id) << " input position: (" << hitPoint.x() / CLHEP::cm
0233                           << ", " << hitPoint.y() / CLHEP::cm << "); position from ID (" << xy.first << ", "
0234                           << xy.second << ") distance " << diff << " Valid " << valid
0235                           << " Rho = " << hitPoint.perp() / CLHEP::cm;
0236   }
0237   return id;
0238 }
0239 
0240 void HGCScintSD::update(const BeginOfJob* job) {
0241   if (hgcons_ != nullptr) {
0242     geom_mode_ = hgcons_->geomMode();
0243     slopeMin_ = hgcons_->minSlope();
0244     levelT1_ = hgcons_->levelTop(0);
0245     levelT2_ = hgcons_->levelTop(1);
0246     firstLayer_ = hgcons_->firstLayer() - 1;
0247 #ifdef EDM_ML_DEBUG
0248     edm::LogVerbatim("HGCSim") << "HGCScintSD::Initialized with mode " << geom_mode_ << " Slope cut " << slopeMin_
0249                                << " top Level " << levelT1_ << ":" << levelT2_ << " FirstLayer " << firstLayer_;
0250 #endif
0251 
0252     numberingScheme_ = std::make_unique<HGCalNumberingScheme>(*hgcons_, mydet_, nameX_, fileName_);
0253   } else {
0254     throw cms::Exception("Unknown", "HGCScintSD") << "Cannot find HGCalDDDConstants for " << nameX_ << "\n";
0255   }
0256 }
0257 
0258 void HGCScintSD::initRun() {}
0259 
0260 bool HGCScintSD::filterHit(CaloG4Hit* aHit, double time) {
0261   return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > eminHit_));
0262 }
0263 
0264 uint32_t HGCScintSD::setDetUnitId(int layer, int module, int cell, int iz, G4ThreeVector& pos) {
0265   uint32_t id = numberingScheme_ ? numberingScheme_->getUnitID(layer, module, cell, iz, pos, weight_) : 0;
0266   return id;
0267 }
0268 
0269 bool HGCScintSD::isItinFidVolume(const G4ThreeVector& pos) {
0270   if (fiducialCut_) {
0271     return (hgcons_->distFromEdgeTrap(pos.x(), pos.y(), pos.z()) > distanceFromEdge_);
0272   } else {
0273     return true;
0274   }
0275 }
[ 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!