Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: HFNoseSD.cc
0003 // Description: Sensitive Detector class for High Granularity Calorimeter
0004 ///////////////////////////////////////////////////////////////////////////////
0005 
0006 #include "DataFormats/Math/interface/FastMath.h"
0007 #include "DataFormats/ForwardDetId/interface/HFNoseDetId.h"
0008 #include "SimG4CMS/Calo/interface/HFNoseSD.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 <iostream>
0022 #include <fstream>
0023 #include <iomanip>
0024 
0025 //#define EDM_ML_DEBUG
0026 
0027 HFNoseSD::HFNoseSD(const std::string& name,
0028                    const HGCalDDDConstants* hgc,
0029                    const SensitiveDetectorCatalog& clg,
0030                    edm::ParameterSet const& p,
0031                    const SimTrackManager* manager)
0032     : CaloSD(name,
0033              clg,
0034              p,
0035              manager,
0036              (float)(p.getParameter<edm::ParameterSet>("HFNoseSD").getParameter<double>("TimeSliceUnit")),
0037              p.getParameter<edm::ParameterSet>("HFNoseSD").getParameter<bool>("IgnoreTrackID")),
0038       hgcons_(hgc),
0039       slopeMin_(0),
0040       levelT1_(99),
0041       levelT2_(99),
0042       useSimWt_(0),
0043       tan30deg_(std::tan(30.0 * CLHEP::deg)) {
0044   numberingScheme_.reset(nullptr);
0045   guardRing_.reset(nullptr);
0046   mouseBite_.reset(nullptr);
0047 
0048   edm::ParameterSet m_HFN = p.getParameter<edm::ParameterSet>("HFNoseSD");
0049   eminHit_ = m_HFN.getParameter<double>("EminHit") * CLHEP::MeV;
0050   fiducialCut_ = m_HFN.getParameter<bool>("FiducialCut");
0051   distanceFromEdge_ = m_HFN.getParameter<double>("DistanceFromEdge");
0052   storeAllG4Hits_ = m_HFN.getParameter<bool>("StoreAllG4Hits");
0053   rejectMB_ = m_HFN.getParameter<bool>("RejectMouseBite");
0054   waferRot_ = m_HFN.getParameter<bool>("RotatedWafer");
0055   cornerMinMask_ = m_HFN.getParameter<int>("CornerMinMask");
0056   angles_ = m_HFN.getUntrackedParameter<std::vector<double>>("WaferAngles");
0057 
0058   nameX_ = ((name.find("HFNoseHits") != std::string::npos) ? "HGCalHFNoseSensitive" : "HFNoseSensitive");
0059 
0060   if (storeAllG4Hits_) {
0061     setUseMap(false);
0062     setNumberCheckedHits(0);
0063   }
0064 
0065 #ifdef EDM_ML_DEBUG
0066   edm::LogVerbatim("HFNSim") << "**************************************************"
0067                              << "\n"
0068                              << "*                                                *"
0069                              << "\n"
0070                              << "* Constructing a HFNoseSD  with name " << name << "\n"
0071                              << "*                                                *"
0072                              << "\n"
0073                              << "**************************************************";
0074 #endif
0075   edm::LogVerbatim("HFNSim") << "HFNoseSD:: Threshold for storing hits: " << eminHit_ << " for " << name;
0076   edm::LogVerbatim("HFNSim") << "Flag for storing individual Geant4 Hits " << storeAllG4Hits_;
0077   edm::LogVerbatim("HFNSim") << "Fiducial volume cut with cut from eta/phi "
0078                              << "boundary " << fiducialCut_ << " at " << distanceFromEdge_;
0079   edm::LogVerbatim("HFNSim") << "Reject MosueBite Flag: " << rejectMB_ << " cuts along " << angles_.size()
0080                              << " axes: " << angles_[0] << ", " << angles_[1];
0081 }
0082 
0083 double HFNoseSD::getEnergyDeposit(const G4Step* aStep) {
0084   double r = aStep->GetPreStepPoint()->GetPosition().perp();
0085   double z = std::abs(aStep->GetPreStepPoint()->GetPosition().z());
0086 #ifdef EDM_ML_DEBUG
0087   G4int parCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
0088   G4String parName = aStep->GetTrack()->GetDefinition()->GetParticleName();
0089   G4LogicalVolume* lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
0090   edm::LogVerbatim("HFNSim") << "HFNoseSD: Hit from standard path from " << lv->GetName() << " for Track "
0091                              << aStep->GetTrack()->GetTrackID() << " (" << parCode << ":" << parName << ") R = " << r
0092                              << " Z = " << z << " slope = " << r / z << ":" << slopeMin_;
0093 #endif
0094   // Apply fiducial cut
0095   if (r < z * slopeMin_) {
0096 #ifdef EDM_ML_DEBUG
0097     edm::LogVerbatim("HFNSim") << "HFNoseSD: Fiducial Volume cut";
0098 #endif
0099     return 0.0;
0100   }
0101 
0102   double wt1 = getResponseWt(aStep->GetTrack());
0103   double wt2 = aStep->GetTrack()->GetWeight();
0104   double destep = weight_ * wt1 * (aStep->GetTotalEnergyDeposit());
0105   if (wt2 > 0)
0106     destep *= wt2;
0107 #ifdef EDM_ML_DEBUG
0108   edm::LogVerbatim("HFNSim") << "HFNoseSD: weights= " << weight_ << ":" << wt1 << ":" << wt2 << " Total weight "
0109                              << weight_ * wt1 * wt2 << " deStep: " << aStep->GetTotalEnergyDeposit() << ":" << destep;
0110 #endif
0111   return destep;
0112 }
0113 
0114 uint32_t HFNoseSD::setDetUnitId(const G4Step* aStep) {
0115   const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0116   const G4VTouchable* touch = preStepPoint->GetTouchable();
0117 
0118   //determine the exact position in global coordinates in the mass geometry
0119   G4ThreeVector hitPoint = preStepPoint->GetPosition();
0120   float globalZ = touch->GetTranslation(0).z();
0121   int iz(globalZ > 0 ? 1 : -1);
0122 
0123   int layer(-1), moduleLev(-1), cell(-1);
0124   if (useSimWt_ > 0) {
0125     layer = touch->GetReplicaNumber(2);
0126     moduleLev = 1;
0127   } else if ((touch->GetHistoryDepth() == levelT1_) || (touch->GetHistoryDepth() == levelT2_)) {
0128     layer = touch->GetReplicaNumber(0);
0129   } else {
0130     layer = touch->GetReplicaNumber(3);
0131     cell = touch->GetReplicaNumber(1);
0132     moduleLev = 2;
0133   }
0134   int module = (moduleLev >= 0) ? touch->GetReplicaNumber(moduleLev) : -1;
0135 #ifdef EDM_ML_DEBUG
0136   edm::LogVerbatim("HFNSim") << "DepthsInside: " << touch->GetHistoryDepth() << " name "
0137                              << touch->GetVolume(0)->GetName() << " layer:module:cell " << layer << ":" << moduleLev
0138                              << ":" << module << ":" << cell;
0139   G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
0140   edm::LogVerbatim("HFNSim") << "Depths: " << touch->GetHistoryDepth() << " name " << touch->GetVolume(0)->GetName()
0141                              << ":" << touch->GetReplicaNumber(0) << "   " << touch->GetVolume(1)->GetName() << ":"
0142                              << touch->GetReplicaNumber(1) << "   " << touch->GetVolume(2)->GetName() << ":"
0143                              << touch->GetReplicaNumber(2) << "   " << touch->GetVolume(3)->GetName() << ":"
0144                              << touch->GetReplicaNumber(3) << "   " << touch->GetVolume(4)->GetName() << ":"
0145                              << touch->GetReplicaNumber(4) << "   "
0146                              << " layer:module:cell " << layer << ":" << module << ":" << cell << " Material "
0147                              << mat->GetName() << ":" << mat->GetRadlen();
0148 #endif
0149   // The following statement should be examined later before elimination
0150   if (aStep->GetPreStepPoint()->GetMaterial()->GetRadlen() > 100000.)
0151     return 0;
0152 
0153   uint32_t id = setDetUnitId(layer, module, cell, iz, hitPoint);
0154   if ((rejectMB_ || fiducialCut_) && id != 0) {
0155     auto uv = HFNoseDetId(id).waferUV();
0156 #ifdef EDM_ML_DEBUG
0157     edm::LogVerbatim("HFNSim") << "ID " << std::hex << id << std::dec << " " << HFNoseDetId(id);
0158 #endif
0159     if ((rejectMB_) && (mouseBite_->exclude(hitPoint, iz, layer, uv.first, uv.second))) {
0160       id = 0;
0161 #ifdef EDM_ML_DEBUG
0162       edm::LogVerbatim("HFNSim") << "Rejected by MouseBite cutoff *****";
0163 #endif
0164     }
0165   }
0166   return id;
0167 }
0168 
0169 void HFNoseSD::update(const BeginOfJob* job) {
0170   if (hgcons_ != nullptr) {
0171     geom_mode_ = hgcons_->geomMode();
0172     slopeMin_ = hgcons_->minSlope();
0173     levelT1_ = hgcons_->levelTop(0);
0174     levelT2_ = hgcons_->levelTop(1);
0175     int useOffset = hgcons_->getParameter()->useOffset_;
0176     double waferSize = hgcons_->waferSize(false);
0177     double mouseBite = hgcons_->mouseBite(false);
0178     mouseBiteCut_ = waferSize * tan30deg_ - mouseBite;
0179     if (useOffset > 0) {
0180       rejectMB_ = true;
0181       fiducialCut_ = true;
0182     }
0183 #ifdef EDM_ML_DEBUG
0184     edm::LogVerbatim("HFNSim") << "HFNoseSD::Initialized with mode " << geom_mode_ << " Slope cut " << slopeMin_
0185                                << " top Level " << levelT1_ << ":" << levelT2_ << " wafer " << waferSize << ":"
0186                                << mouseBite << " useOffset " << useOffset;
0187 #endif
0188 
0189     numberingScheme_ = std::make_unique<HFNoseNumberingScheme>(*hgcons_);
0190     if (rejectMB_)
0191       mouseBite_ = std::make_unique<HGCMouseBite>(*hgcons_, angles_, mouseBiteCut_, waferRot_);
0192   } else {
0193     throw cms::Exception("Unknown", "HFNoseSD") << "Cannot find HGCalDDDConstants for " << nameX_ << "\n";
0194   }
0195 }
0196 
0197 void HFNoseSD::initRun() {}
0198 
0199 bool HFNoseSD::filterHit(CaloG4Hit* aHit, double time) {
0200   return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > eminHit_));
0201 }
0202 
0203 uint32_t HFNoseSD::setDetUnitId(int layer, int module, int cell, int iz, G4ThreeVector& pos) {
0204   uint32_t id = numberingScheme_ ? numberingScheme_->getUnitID(layer, module, cell, iz, pos, weight_) : 0;
0205   if (cornerMinMask_ > 2) {
0206     if (hgcons_->maskCell(DetId(id), cornerMinMask_))
0207       id = 0;
0208   }
0209   return id;
0210 }
0211 
0212 bool HFNoseSD::isItinFidVolume(const G4ThreeVector& pos) {
0213   if (fiducialCut_) {
0214     return (hgcons_->distFromEdgeHex(pos.x(), pos.y(), pos.z()) > distanceFromEdge_);
0215   } else {
0216     return true;
0217   }
0218 }
[ 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!