Calo/src/HGCScintSD.cc

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 "SimG4Core/Notification/interface/TrackInformation.h"
0011 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0012 #include "FWCore/Utilities/interface/Exception.h"
0013 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0014 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
0015 #include "G4LogicalVolumeStore.hh"
0016 #include "G4LogicalVolume.hh"
0017 #include "G4Step.hh"
0018 #include "G4Track.hh"
0019 #include "G4ParticleTable.hh"
0020 #include "G4VProcess.hh"
0021 #include "G4Trap.hh"
0022
0023 #include <fstream>
0024 #include <iomanip>
0025 #include <iostream>
0026 #include <memory>
0027
0028 //#define EDM_ML_DEBUG
0029
0030 HGCScintSD::HGCScintSD(const std::string& name,
0031                        const HGCalDDDConstants* hgc,
0032                        const SensitiveDetectorCatalog& clg,
0033                        edm::ParameterSet const& p,
0034                        const SimTrackManager* manager)
0035     : CaloSD(name,
0036              clg,
0037              p,
0038              manager,
0039              (float)(p.getParameter<edm::ParameterSet>("HGCSD").getParameter<double>("TimeSliceUnit")),
0040              p.getParameter<edm::ParameterSet>("HGCSD").getParameter<bool>("IgnoreTrackID")),
0041       hgcons_(hgc),
0042       slopeMin_(0),
0043       levelT1_(99),
0044       levelT2_(99) {
0045   numberingScheme_.reset(nullptr);
0046
0047   edm::ParameterSet m_HGC = p.getParameter<edm::ParameterSet>("HGCScintSD");
0048   eminHit_ = m_HGC.getParameter<double>("EminHit") * CLHEP::MeV;
0049   fiducialCut_ = m_HGC.getParameter<bool>("FiducialCut");
0050   distanceFromEdge_ = m_HGC.getParameter<double>("DistanceFromEdge");
0051   useBirk_ = m_HGC.getParameter<bool>("UseBirkLaw");
0052   birk1_ = m_HGC.getParameter<double>("BirkC1") * (CLHEP::g / (CLHEP::MeV * CLHEP::cm2));
0053   birk2_ = m_HGC.getParameter<double>("BirkC2");
0054   birk3_ = m_HGC.getParameter<double>("BirkC3");
0055   storeAllG4Hits_ = m_HGC.getParameter<bool>("StoreAllG4Hits");
0056
0057   if (storeAllG4Hits_) {
0058     setUseMap(false);
0059     setNumberCheckedHits(0);
0060   }
0061
0062   //this is defined in the hgcsens.xml
0063   G4String myName = name;
0064   mydet_ = DetId::Forward;
0065   nameX_ = "HGCal";
0066   if (myName.find("HitsHEback") != std::string::npos) {
0067     mydet_ = DetId::HGCalHSc;
0068     nameX_ = "HGCalHEScintillatorSensitive";
0069   }
0070
0071 #ifdef EDM_ML_DEBUG
0072   edm::LogVerbatim("HGCSim") << "**************************************************"
0073                              << "\n"
0074                              << "*                                                *"
0075                              << "\n"
0076                              << "* Constructing a HGCScintSD  with name " << name << "\n"
0077                              << "*                                                *"
0078                              << "\n"
0079                              << "**************************************************";
0080 #endif
0081   edm::LogVerbatim("HGCSim") << "HGCScintSD:: Threshold for storing hits: " << eminHit_ << " for " << nameX_
0082                              << " detector " << mydet_;
0083   edm::LogVerbatim("HGCSim") << "Flag for storing individual Geant4 Hits " << storeAllG4Hits_;
0084   edm::LogVerbatim("HGCSim") << "Fiducial volume cut with cut from eta/phi "
0085                              << "boundary " << fiducialCut_ << " at " << distanceFromEdge_;
0086   edm::LogVerbatim("HGCSim") << "Use of Birks law is set to      " << useBirk_
0087                              << "  with three constants kB = " << birk1_ << ", C1 = " << birk2_ << ", C2 = " << birk3_;
0088 }
0089
0090 double HGCScintSD::getEnergyDeposit(const G4Step* aStep) {
0091   double r = aStep->GetPreStepPoint()->GetPosition().perp();
0092   double z = std::abs(aStep->GetPreStepPoint()->GetPosition().z());
0093 #ifdef EDM_ML_DEBUG
0094   G4int parCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
0095   G4String parName = aStep->GetTrack()->GetDefinition()->GetParticleName();
0096   G4LogicalVolume* lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
0097   edm::LogVerbatim("HGCSim") << "HGCScintSD: Hit from standard path from " << lv->GetName() << " for Track "
0098                              << aStep->GetTrack()->GetTrackID() << " (" << parCode << ":" << parName << ") R = " << r
0099                              << " Z = " << z << " slope = " << r / z << ":" << slopeMin_;
0100 #endif
0101   // Apply fiducial cut
0102   if (r < z * slopeMin_) {
0103 #ifdef EDM_ML_DEBUG
0104     edm::LogVerbatim("HGCSim") << "HGCScintSD: Fiducial Volume cut";
0105 #endif
0106     return 0.0;
0107   }
0108
0109   double wt1 = getResponseWt(aStep->GetTrack());
0110   double wt2 = aStep->GetTrack()->GetWeight();
0111   double wt3 = (useBirk_ ? getAttenuation(aStep, birk1_, birk2_, birk3_) : 1.0);
0112   double destep = weight_ * wt1 * wt3 * (aStep->GetTotalEnergyDeposit());
0113   if (wt2 > 0)
0114     destep *= wt2;
0115 #ifdef EDM_ML_DEBUG
0116   edm::LogVerbatim("HGCalSim") << "HGCScintSD: weights= " << weight_ << ":" << wt1 << ":" << wt2 << ":" << wt3
0117                                << " Total weight " << weight_ * wt1 * wt2 * wt3
0118                                << " deStep: " << aStep->GetTotalEnergyDeposit() << ":" << destep;
0119 #endif
0120   return destep;
0121 }
0122
0123 uint32_t HGCScintSD::setDetUnitId(const G4Step* aStep) {
0124   const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0125   const G4VTouchable* touch = preStepPoint->GetTouchable();
0126
0127 #ifdef EDM_ML_DEBUG
0128   edm::LogVerbatim("HGCSim") << "DepthsTop: " << touch->GetHistoryDepth() << ":" << levelT1_ << ":" << levelT2_;
0129   printDetectorLevels(touch);
0130 #endif
0131   //determine the exact position in global coordinates in the mass geometry
0132   G4ThreeVector hitPoint = preStepPoint->GetPosition();
0133   float globalZ = touch->GetTranslation(0).z();
0134   int iz(globalZ > 0 ? 1 : -1);
0135
0136   int layer(0), module(-1), cell(-1);
0137   if ((geom_mode_ == HGCalGeometryMode::TrapezoidModule) || (geom_mode_ == HGCalGeometryMode::TrapezoidCassette)) {
0138     layer = touch->GetReplicaNumber(1);
0139   } else if ((touch->GetHistoryDepth() == levelT1_) || (touch->GetHistoryDepth() == levelT2_)) {
0140     layer = touch->GetReplicaNumber(0);
0141 #ifdef EDM_ML_DEBUG
0142     edm::LogVerbatim("HGCSim") << "DepthsTop: " << touch->GetHistoryDepth() << ":" << levelT1_ << ":" << levelT2_
0143                                << " name " << touch->GetVolume(0)->GetName() << " layer:module:cell " << layer << ":"
0144                                << module << ":" << cell;
0145 #endif
0146   } else {
0147     layer = touch->GetReplicaNumber(3);
0148     module = touch->GetReplicaNumber(2);
0149     cell = touch->GetReplicaNumber(1);
0150 #ifdef EDM_ML_DEBUG
0151     edm::LogVerbatim("HGCSim") << "DepthsInside: " << touch->GetHistoryDepth() << " name "
0152                                << touch->GetVolume(0)->GetName() << " layer:module:cell " << layer << ":" << module
0153                                << ":" << cell;
0154 #endif
0155   }
0156 #ifdef EDM_ML_DEBUG
0157   G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
0158   edm::LogVerbatim("HGCSim") << "Depths: " << touch->GetHistoryDepth() << " name " << touch->GetVolume(0)->GetName()
0159                              << ":" << touch->GetReplicaNumber(0) << "   " << touch->GetVolume(1)->GetName() << ":"
0160                              << touch->GetReplicaNumber(1) << "   " << touch->GetVolume(2)->GetName() << ":"
0161                              << touch->GetReplicaNumber(2) << "   " << touch->GetVolume(3)->GetName() << ":"
0162                              << touch->GetReplicaNumber(3) << "   " << touch->GetVolume(4)->GetName() << ":"
0163                              << touch->GetReplicaNumber(4) << "   "
0164                              << " layer:module:cell " << layer << ":" << module << ":" << cell << " Material "
0165                              << mat->GetName() << ":" << mat->GetRadlen();
0166 #endif
0167   // The following statement should be examined later before elimination
0168   if (aStep->GetPreStepPoint()->GetMaterial()->GetRadlen() > 100000.)
0169     return 0;
0170
0171   uint32_t id = setDetUnitId(layer, module, cell, iz, hitPoint);
0172   if (!isItinFidVolume(hitPoint)) {
0173 #ifdef EDM_ML_DEBUG
0174     edm::LogVerbatim("HGCSim") << "ID " << std::hex << id << std::dec << " " << HGCScintillatorDetId(id)
0175                                << " is rejected by fiducilal volume cut";
0176 #endif
0177     id = 0;
0178   }
0179   return id;
0180 }
0181
0182 void HGCScintSD::update(const BeginOfJob* job) {
0183   if (hgcons_ != nullptr) {
0184     geom_mode_ = hgcons_->geomMode();
0185     slopeMin_ = hgcons_->minSlope();
0186     levelT1_ = hgcons_->levelTop(0);
0187     levelT2_ = hgcons_->levelTop(1);
0188 #ifdef EDM_ML_DEBUG
0189     edm::LogVerbatim("HGCSim") << "HGCScintSD::Initialized with mode " << geom_mode_ << " Slope cut " << slopeMin_
0190                                << " top Level " << levelT1_ << ":" << levelT2_;
0191 #endif
0192
0193     numberingScheme_ = std::make_unique<HGCalNumberingScheme>(*hgcons_, mydet_, nameX_);
0194   } else {
0195     throw cms::Exception("Unknown", "HGCScintSD") << "Cannot find HGCalDDDConstants for " << nameX_ << "\n";
0196   }
0197 }
0198
0199 void HGCScintSD::initRun() {}
0200
0201 bool HGCScintSD::filterHit(CaloG4Hit* aHit, double time) {
0202   return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > eminHit_));
0203 }
0204
0205 uint32_t HGCScintSD::setDetUnitId(int layer, int module, int cell, int iz, G4ThreeVector& pos) {
0206   uint32_t id = numberingScheme_ ? numberingScheme_->getUnitID(layer, module, cell, iz, pos, weight_) : 0;
0207   return id;
0208 }
0209
0210 bool HGCScintSD::isItinFidVolume(const G4ThreeVector& pos) {
0211   if (fiducialCut_) {
0212     return (hgcons_->distFromEdgeTrap(pos.x(), pos.y(), pos.z()) > distanceFromEdge_);
0213   } else {
0214     return true;
0215   }
0216 }