Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 10:03:59 

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: ECalSD.cc
0003 // Description: Sensitive Detector class for electromagnetic calorimeters
0004 ///////////////////////////////////////////////////////////////////////////////
0005 #include "SimG4CMS/Calo/interface/ECalSD.h"
0006 #include "SimG4CMS/Calo/interface/EcalDumpGeometry.h"
0007 #include "SimG4Core/Notification/interface/TrackInformation.h"
0008 #include "Geometry/EcalCommonData/interface/EcalBarrelNumberingScheme.h"
0009 #include "Geometry/EcalCommonData/interface/EcalBaseNumber.h"
0010 #include "Geometry/EcalCommonData/interface/EcalEndcapNumberingScheme.h"
0011 #include "Geometry/EcalCommonData/interface/EcalPreshowerNumberingScheme.h"
0012 #include "Geometry/EcalCommonData/interface/ESTBNumberingScheme.h"
0013 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0014 #include "DataFormats/Math/interface/GeantUnits.h"
0015 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0016 
0017 #include "FWCore/Framework/interface/ESTransientHandle.h"
0018 #include "FWCore/ServiceRegistry/interface/Service.h"
0019 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0020 
0021 #include "DD4hep/Filter.h"
0022 
0023 #include "G4LogicalVolumeStore.hh"
0024 #include "G4LogicalVolume.hh"
0025 #include "G4Step.hh"
0026 #include "G4Track.hh"
0027 #include "G4VProcess.hh"
0028 #include "G4SystemOfUnits.hh"
0029 
0030 #include <algorithm>
0031 
0032 //#define EDM_ML_DEBUG
0033 
0034 using namespace geant_units::operators;
0035 
0036 template <class T>
0037 bool any(const std::vector<T>& v, const T& what) {
0038   return std::find(v.begin(), v.end(), what) != v.end();
0039 }
0040 
0041 ECalSD::ECalSD(const std::string& name,
0042                const EcalSimulationParameters* ecpar,
0043                const SensitiveDetectorCatalog& clg,
0044                edm::ParameterSet const& p,
0045                const SimTrackManager* manager)
0046     : CaloSD(name,
0047              clg,
0048              p,
0049              manager,
0050              (float)(p.getParameter<edm::ParameterSet>("ECalSD").getParameter<double>("TimeSliceUnit")),
0051              p.getParameter<edm::ParameterSet>("ECalSD").getParameter<bool>("IgnoreTrackID")),
0052       ecalSimParameters_(ecpar) {
0053   numberingScheme_.reset(nullptr);
0054   //   static SimpleConfigurable<bool>   on1(false,  "ECalSD:UseBirkLaw");
0055   //   static SimpleConfigurable<double> bk1(0.00463,"ECalSD:BirkC1");
0056   //   static SimpleConfigurable<double> bk2(-0.03,  "ECalSD:BirkC2");
0057   //   static SimpleConfigurable<double> bk3(1.0,    "ECalSD:BirkC3");
0058   // Values from NIM A484 (2002) 239-244: as implemented in Geant3
0059   //   useBirk          = on1.value();
0060   //   birk1            = bk1.value()*(g/(MeV*cm2));
0061   //   birk2            = bk2.value()*(g/(MeV*cm2))*(g/(MeV*cm2));
0062   edm::ParameterSet m_EC = p.getParameter<edm::ParameterSet>("ECalSD");
0063   useBirk = m_EC.getParameter<bool>("UseBirkLaw");
0064   useBirkL3 = m_EC.getParameter<bool>("BirkL3Parametrization");
0065   double bunit = (CLHEP::g / (CLHEP::MeV * CLHEP::cm2));
0066   birk1 = m_EC.getParameter<double>("BirkC1") * bunit;
0067   birk2 = m_EC.getParameter<double>("BirkC2");
0068   birk3 = m_EC.getParameter<double>("BirkC3");
0069   birkSlope = m_EC.getParameter<double>("BirkSlope");
0070   birkCut = m_EC.getParameter<double>("BirkCut");
0071   slopeLY = m_EC.getParameter<double>("SlopeLightYield");
0072   storeTrack = m_EC.getParameter<bool>("StoreSecondary");
0073   crystalMat = m_EC.getUntrackedParameter<std::string>("XtalMat", "E_PbWO4");
0074   bool isItTB = m_EC.getUntrackedParameter<bool>("TestBeam", false);
0075   bool nullNS = m_EC.getUntrackedParameter<bool>("NullNumbering", false);
0076   storeRL = m_EC.getUntrackedParameter<bool>("StoreRadLength", false);
0077   scaleRL = m_EC.getUntrackedParameter<double>("ScaleRadLength", 1.0);
0078   int dumpGeom = m_EC.getUntrackedParameter<int>("DumpGeometry", 0);
0079 
0080   //Changes for improved timing simulation
0081   storeLayerTimeSim = m_EC.getUntrackedParameter<bool>("StoreLayerTimeSim", false);
0082 
0083   ageingWithSlopeLY = m_EC.getUntrackedParameter<bool>("AgeingWithSlopeLY", false);
0084   if (ageingWithSlopeLY)
0085     ageing.setLumies(p.getParameter<edm::ParameterSet>("ECalSD").getParameter<double>("DelivLuminosity"),
0086                      p.getParameter<edm::ParameterSet>("ECalSD").getParameter<double>("InstLuminosity"));
0087 
0088   if (ecalSimParameters_ == nullptr) {
0089     edm::LogError("EcalSim") << "ECalSD : Cannot find EcalSimulationParameters for " << name;
0090     throw cms::Exception("Unknown", "ECalSD") << "Cannot find EcalSimulationParameters for " << name << "\n";
0091   }
0092 
0093   // Use of Weight
0094   useWeight = ecalSimParameters_->useWeight_;
0095 #ifdef EDM_ML_DEBUG
0096   edm::LogVerbatim("EcalSim") << "ECalSD:: useWeight " << useWeight;
0097 #endif
0098   depth1Name = ecalSimParameters_->depth1Name_;
0099   depth2Name = ecalSimParameters_->depth2Name_;
0100 #ifdef EDM_ML_DEBUG
0101   edm::LogVerbatim("EcalSim") << "Names (Depth 1):" << depth1Name << " (Depth 2):" << depth2Name << std::endl;
0102 #endif
0103   int type(-1);
0104   bool dump(false);
0105   EcalNumberingScheme* scheme = nullptr;
0106   if (nullNS) {
0107     scheme = nullptr;
0108   } else if (name == "EcalHitsEB") {
0109     scheme = dynamic_cast<EcalNumberingScheme*>(new EcalBarrelNumberingScheme());
0110     type = 0;
0111     dump = ((dumpGeom % 10) > 0);
0112   } else if (name == "EcalHitsEE") {
0113     scheme = dynamic_cast<EcalNumberingScheme*>(new EcalEndcapNumberingScheme());
0114     type = 1;
0115     dump = (((dumpGeom / 10) % 10) > 0);
0116   } else if (name == "EcalHitsES") {
0117     if (isItTB)
0118       scheme = dynamic_cast<EcalNumberingScheme*>(new ESTBNumberingScheme());
0119     else
0120       scheme = dynamic_cast<EcalNumberingScheme*>(new EcalPreshowerNumberingScheme());
0121     useWeight = false;
0122     type = 2;
0123     dump = (((dumpGeom / 100) % 10) > 0);
0124   } else {
0125     edm::LogWarning("EcalSim") << "ECalSD: ReadoutName not supported";
0126   }
0127   int type0 = dumpGeom / 1000;
0128   type += (10 * type0);
0129 
0130   if (nullptr != scheme)
0131     setNumberingScheme(scheme);
0132 #ifdef EDM_ML_DEBUG
0133   edm::LogVerbatim("EcalSim") << "Constructing a ECalSD  with name " << GetName();
0134 #endif
0135   if (useWeight) {
0136     edm::LogVerbatim("EcalSim") << "ECalSD:: Use of Birks law is set to " << useBirk
0137                                 << " with three constants kB = " << birk1 / bunit << ", C1 = " << birk2
0138                                 << ", C2 = " << birk3 << "\n         Use of L3 parametrization " << useBirkL3
0139                                 << " with slope " << birkSlope << " and cut off " << birkCut << "\n"
0140                                 << "         Slope for Light yield is set to " << slopeLY;
0141   } else {
0142     edm::LogVerbatim("EcalSim") << "ECalSD:: energy deposit is not corrected "
0143                                 << " by Birk or light yield curve";
0144   }
0145 
0146   edm::LogVerbatim("EcalSim") << "ECalSD:: Suppression Flag " << suppressHeavy << "\tprotons below "
0147                               << kmaxProton / CLHEP::MeV << " MeV,\tneutrons below " << kmaxNeutron / CLHEP::MeV
0148                               << " MeV,\tions below " << kmaxIon / CLHEP::MeV << " MeV \n\tDepth1 Name = " << depth1Name
0149                               << "\tDepth2 Name = " << depth2Name << "\n\tstoreRL " << storeRL << ":" << scaleRL
0150                               << "\tstoreLayerTimeSim " << storeLayerTimeSim << "\n\ttime Granularity "
0151                               << p.getParameter<edm::ParameterSet>("ECalSD").getParameter<double>("TimeSliceUnit")
0152                               << " ns";
0153   if (useWeight)
0154     initMap();
0155 #ifdef plotDebug
0156   edm::Service<TFileService> tfile;
0157   if (tfile.isAvailable()) {
0158     TFileDirectory ecDir = tfile->mkdir("ProfileFromECalSD");
0159     static const std::string ctype[4] = {"EB", "EBref", "EE", "EERef"};
0160     for (int k = 0; k < 4; ++k) {
0161       std::string name = "ECLL_" + ctype[k];
0162       std::string title = "Local vs Global for " + ctype[k];
0163       double xmin = (k > 1) ? 3000.0 : 1000.0;
0164       g2L_[k] = ecDir.make<TH2F>(name.c_str(), title.c_str(), 100, xmin, xmin + 1000., 100, 0.0, 3000.);
0165     }
0166   } else {
0167     for (int k = 0; k < 4; ++k)
0168       g2L_[k] = 0;
0169   }
0170 #endif
0171   if (dump) {
0172     const auto& lvNames = clg.logicalNames(name);
0173     EcalDumpGeometry geom(lvNames, depth1Name, depth2Name, type);
0174     geom.update();
0175   }
0176 }
0177 
0178 ECalSD::~ECalSD() {}
0179 
0180 double ECalSD::getEnergyDeposit(const G4Step* aStep) {
0181   const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0182   const G4Track* theTrack = aStep->GetTrack();
0183   double edep = aStep->GetTotalEnergyDeposit();
0184 
0185   // take into account light collection curve for crystals
0186   double weight = 1.;
0187   if (suppressHeavy) {
0188     TrackInformation* trkInfo = (TrackInformation*)(theTrack->GetUserInformation());
0189     if (trkInfo) {
0190       int pdg = theTrack->GetDefinition()->GetPDGEncoding();
0191       if (!(trkInfo->isPrimary())) {  // Only secondary particles
0192         double ke = theTrack->GetKineticEnergy();
0193         if (((pdg / 1000000000 == 1 && ((pdg / 10000) % 100) > 0 && ((pdg / 10) % 100) > 0)) && (ke < kmaxIon))
0194           weight = 0;
0195         if ((pdg == 2212) && (ke < kmaxProton))
0196           weight = 0;
0197         if ((pdg == 2112) && (ke < kmaxNeutron))
0198           weight = 0;
0199       }
0200     }
0201   }
0202   const G4LogicalVolume* lv = preStepPoint->GetTouchable()->GetVolume(0)->GetLogicalVolume();
0203   double wt1 = 1.0;
0204   if (useWeight && !any(noWeight, lv)) {
0205     weight *= curve_LY(lv);
0206     if (useBirk) {
0207       if (useBirkL3)
0208         weight *= getBirkL3(aStep);
0209       else
0210         weight *= getAttenuation(aStep, birk1, birk2, birk3);
0211     }
0212     wt1 = getResponseWt(theTrack);
0213   }
0214   edep *= weight * wt1;
0215   // Russian Roulette
0216   double wt2 = theTrack->GetWeight();
0217   if (wt2 > 0.0) {
0218     edep *= wt2;
0219   }
0220 #ifdef EDM_ML_DEBUG
0221   edm::LogVerbatim("EcalSim") << lv->GetName() << " " << dd4hep::dd::noNamespace(lv->GetName())
0222                               << " Light Collection Efficiency " << weight << ":" << wt1 << " wt2= " << wt2
0223                               << " Weighted Energy Deposit " << edep / CLHEP::MeV << " MeV at "
0224                               << preStepPoint->GetPosition();
0225 #endif
0226   return edep;
0227 }
0228 
0229 double ECalSD::EnergyCorrected(const G4Step& step, const G4Track* track) {
0230   const G4StepPoint* hitPoint = step.GetPreStepPoint();
0231   const G4LogicalVolume* lv = hitPoint->GetTouchable()->GetVolume(0)->GetLogicalVolume();
0232   if (lv->GetSensitiveDetector() != this)
0233     return 0.0;
0234 
0235   double edep = step.GetTotalEnergyDeposit();
0236   if (useWeight && !any(noWeight, lv)) {
0237     currentLocalPoint = setToLocal(hitPoint->GetPosition(), hitPoint->GetTouchable());
0238     auto ite = xtalLMap.find(lv);
0239     crystalLength = (ite == xtalLMap.end()) ? 230._mm : std::abs(ite->second);
0240     crystalDepth = (ite == xtalLMap.end()) ? 0.0 : (std::abs(0.5 * (ite->second) + currentLocalPoint.z()));
0241     edep *= curve_LY(lv) * getResponseWt(track);
0242   }
0243   return edep;
0244 }
0245 
0246 int ECalSD::getTrackID(const G4Track* aTrack) {
0247   int primaryID(0);
0248   if (storeTrack && depth > 0) {
0249     forceSave = true;
0250     primaryID = aTrack->GetTrackID();
0251   } else {
0252     primaryID = CaloSD::getTrackID(aTrack);
0253   }
0254   return primaryID;
0255 }
0256 
0257 uint16_t ECalSD::getDepth(const G4Step* aStep) {
0258   // this method should be called first at a step
0259   const G4StepPoint* hitPoint = aStep->GetPreStepPoint();
0260   currentLocalPoint = setToLocal(hitPoint->GetPosition(), hitPoint->GetTouchable());
0261   const G4LogicalVolume* lv = hitPoint->GetTouchable()->GetVolume(0)->GetLogicalVolume();
0262 
0263   auto ite = xtalLMap.find(lv);
0264   crystalLength = (ite == xtalLMap.end()) ? 230._mm : std::abs(ite->second);
0265   crystalDepth = (ite == xtalLMap.end()) ? 0.0 : (std::abs(0.5 * (ite->second) + currentLocalPoint.z()));
0266   depth = any(useDepth1, lv) ? 1 : (any(useDepth2, lv) ? 2 : 0);
0267   uint16_t depth1(0), depth2(0);
0268   if (storeRL) {
0269     depth1 = (ite == xtalLMap.end()) ? 0 : (((ite->second) >= 0) ? 0 : PCaloHit::kEcalDepthRefz);
0270     depth2 = getRadiationLength(hitPoint, lv);
0271     depth |= (((depth2 & PCaloHit::kEcalDepthMask) << PCaloHit::kEcalDepthOffset) | depth1);
0272   } else if (storeLayerTimeSim) {
0273     depth2 = getLayerIDForTimeSim();
0274     depth |= ((depth2 & PCaloHit::kEcalDepthMask) << PCaloHit::kEcalDepthOffset);
0275   }
0276 #ifdef EDM_ML_DEBUG
0277   if (isXtal(lv))
0278     edm::LogVerbatim("EcalSimX") << "ECalSD::Volume " << lv->GetName() << " DetId " << std::hex << setDetUnitId(aStep)
0279                                  << std::dec << " Global " << (hitPoint->GetPosition()).rho() << ":"
0280                                  << (hitPoint->GetPosition()).z() << " Local Z " << currentLocalPoint.z() << " Depth "
0281                                  << crystalDepth;
0282   edm::LogVerbatim("EcalSim") << "ECalSD::Depth " << std::hex << depth1 << ":" << depth2 << ":" << depth << std::dec
0283                               << " L " << (ite == xtalLMap.end()) << ":" << ite->second << " local "
0284                               << currentLocalPoint << " Crystal length " << crystalLength << ":" << crystalDepth;
0285 #endif
0286   return depth;
0287 }
0288 
0289 uint16_t ECalSD::getRadiationLength(const G4StepPoint* hitPoint, const G4LogicalVolume* lv) {
0290   uint16_t thisX0 = 0;
0291   if (useWeight) {
0292     double radl = hitPoint->GetMaterial()->GetRadlen();
0293     thisX0 = (uint16_t)floor(scaleRL * crystalDepth / radl);
0294 #ifdef plotDebug
0295     const std::string& lvname = dd4hep::dd::noNamespace(lv->GetName());
0296     int k1 = (lvname.find("EFRY") != std::string::npos) ? 2 : 0;
0297     int k2 = (lvname.find("refl") != std::string::npos) ? 1 : 0;
0298     int kk = k1 + k2;
0299     double rz = (k1 == 0) ? (hitPoint->GetPosition()).rho() : std::abs((hitPoint->GetPosition()).z());
0300     edm::LogVerbatim("EcalSim") << lvname << " # " << k1 << ":" << k2 << ":" << kk << " rz " << rz << " D " << thisX0;
0301     g2L_[kk]->Fill(rz, thisX0);
0302 #endif
0303 #ifdef EDM_ML_DEBUG
0304     G4ThreeVector localPoint = setToLocal(hitPoint->GetPosition(), hitPoint->GetTouchable());
0305     edm::LogVerbatim("EcalSim") << lv->GetName() << " " << dd4hep::dd::noNamespace(lv->GetName()) << " Global "
0306                                 << hitPoint->GetPosition() << ":" << (hitPoint->GetPosition()).rho() << " Local "
0307                                 << localPoint << " Crystal Length " << crystalLength << " Radl " << radl
0308                                 << " crystalDepth " << crystalDepth << " Index " << thisX0 << " : "
0309                                 << getLayerIDForTimeSim();
0310 #endif
0311   }
0312   return thisX0;
0313 }
0314 
0315 uint16_t ECalSD::getLayerIDForTimeSim() {
0316   const double invLayerSize = 0.1;  //layer size in 1/mm
0317   return (int)crystalDepth * invLayerSize;
0318 }
0319 
0320 uint32_t ECalSD::setDetUnitId(const G4Step* aStep) {
0321   if (numberingScheme_.get() == nullptr) {
0322     return EBDetId(1, 1)();
0323   } else {
0324     getBaseNumber(aStep);
0325     return numberingScheme_->getUnitID(theBaseNumber);
0326   }
0327 }
0328 
0329 void ECalSD::setNumberingScheme(EcalNumberingScheme* scheme) {
0330   if (scheme != nullptr) {
0331     edm::LogVerbatim("EcalSim") << "EcalSD: updates numbering scheme for " << GetName();
0332     numberingScheme_.reset(scheme);
0333   }
0334 }
0335 
0336 void ECalSD::initMap() {
0337   std::vector<const G4LogicalVolume*> lvused;
0338   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0339   std::map<const std::string, const G4LogicalVolume*> nameMap;
0340   for (auto lvi = lvs->begin(), lve = lvs->end(); lvi != lve; ++lvi)
0341     nameMap.emplace(dd4hep::dd::noNamespace((*lvi)->GetName()), *lvi);
0342 
0343   for (unsigned int it = 0; it < ecalSimParameters_->lvNames_.size(); ++it) {
0344     const std::string& matname = ecalSimParameters_->matNames_[it];
0345     const std::string& lvname = ecalSimParameters_->lvNames_[it];
0346     const G4LogicalVolume* lv = nameMap[lvname];
0347     int ibec = (lvname.find("EFRY") == std::string::npos) ? 0 : 1;
0348     int iref = (lvname.find("refl") == std::string::npos) ? 0 : 1;
0349     int type = (ibec + iref == 1) ? 1 : -1;
0350     if (depth1Name != " ") {
0351       if (strncmp(lvname.c_str(), depth1Name.c_str(), 4) == 0) {
0352         if (!any(useDepth1, lv)) {
0353           useDepth1.push_back(lv);
0354 #ifdef EDM_ML_DEBUG
0355           edm::LogVerbatim("EcalSim") << "ECalSD::initMap Logical Volume " << lvname << " in Depth 1 volume list";
0356 #endif
0357         }
0358         const G4LogicalVolume* lvr = nameMap[lvname + "_refl"];
0359         if (lvr != nullptr && !any(useDepth1, lvr)) {
0360           useDepth1.push_back(lvr);
0361 #ifdef EDM_ML_DEBUG
0362           edm::LogVerbatim("EcalSim") << "ECalSD::initMap Logical Volume " << lvname << "_refl"
0363                                       << " in Depth 1 volume list";
0364 #endif
0365         }
0366       }
0367     }
0368     if (depth2Name != " ") {
0369       if (strncmp(lvname.c_str(), depth2Name.c_str(), 4) == 0) {
0370         if (!any(useDepth2, lv)) {
0371           useDepth2.push_back(lv);
0372 #ifdef EDM_ML_DEBUG
0373           edm::LogVerbatim("EcalSim") << "ECalSD::initMap Logical Volume " << lvname << " in Depth 2 volume list";
0374 #endif
0375         }
0376         const G4LogicalVolume* lvr = nameMap[lvname + "_refl"];
0377         if (lvr != nullptr && !any(useDepth2, lvr)) {
0378           useDepth2.push_back(lvr);
0379 #ifdef EDM_ML_DEBUG
0380           edm::LogVerbatim("EcalSim") << "ECalSD::initMap Logical Volume " << lvname << "_refl"
0381                                       << " in Depth 2 volume list";
0382 #endif
0383         }
0384       }
0385     }
0386     if (lv != nullptr) {
0387       if (crystalMat.size() == matname.size() && !strcmp(crystalMat.c_str(), matname.c_str())) {
0388         if (!any(lvused, lv)) {
0389           lvused.push_back(lv);
0390           double dz = ecalSimParameters_->dzs_[it];
0391           xtalLMap.insert(std::pair<const G4LogicalVolume*, double>(lv, dz * type));
0392           lv = nameMap[lvname + "_refl"];
0393           if (lv != nullptr) {
0394             xtalLMap.insert(std::pair<const G4LogicalVolume*, double>(lv, -dz * type));
0395           }
0396         }
0397       } else {
0398         if (!any(noWeight, lv)) {
0399           noWeight.push_back(lv);
0400 #ifdef EDM_ML_DEBUG
0401           edm::LogVerbatim("EcalSim") << "ECalSD::initMap Logical Volume " << lvname << " Material " << matname
0402                                       << " in noWeight list";
0403 #endif
0404         }
0405         lv = nameMap[lvname];
0406         if (lv != nullptr && !any(noWeight, lv)) {
0407           noWeight.push_back(lv);
0408 #ifdef EDM_ML_DEBUG
0409           edm::LogVerbatim("EcalSim") << "ECalSD::initMap Logical Volume " << lvname << " Material " << matname
0410                                       << " in noWeight list";
0411 #endif
0412         }
0413       }
0414     }
0415   }
0416 #ifdef EDM_ML_DEBUG
0417   edm::LogVerbatim("EcalSim") << "ECalSD: Length Table:";
0418   int i = 0;
0419   for (auto ite : xtalLMap) {
0420     std::string name("Unknown");
0421     if (ite.first != nullptr)
0422       name = dd4hep::dd::noNamespace((ite.first)->GetName());
0423     edm::LogVerbatim("EcalSim") << " " << i << " " << ite.first << " " << name << " L = " << ite.second;
0424     ++i;
0425   }
0426 #endif
0427 }
0428 
0429 double ECalSD::curve_LY(const G4LogicalVolume* lv) {
0430   double weight = 1.;
0431   if (ageingWithSlopeLY) {
0432     //position along the crystal in mm from 0 to 230 (in EB)
0433     if (crystalDepth >= -0.1 || crystalDepth <= crystalLength + 0.1)
0434       weight = ageing.calcLightCollectionEfficiencyWeighted(currentID[0].unitID(), crystalDepth / crystalLength);
0435   } else {
0436     double dapd = crystalLength - crystalDepth;
0437     if (dapd >= -0.1 || dapd <= crystalLength + 0.1) {
0438       if (dapd <= 100.)
0439         weight = 1.0 + slopeLY - dapd * 0.01 * slopeLY;
0440     } else {
0441       edm::LogWarning("EcalSim") << "ECalSD: light coll curve : wrong distance "
0442                                  << "to APD " << dapd << " crlength = " << crystalLength << ":" << crystalDepth
0443                                  << " crystal name = " << lv->GetName() << " " << dd4hep::dd::noNamespace(lv->GetName())
0444                                  << " z of localPoint = " << currentLocalPoint.z() << " take weight = " << weight;
0445     }
0446   }
0447 #ifdef EDM_ML_DEBUG
0448   edm::LogVerbatim("EcalSim") << "ECalSD: light coll curve : crlength = " << crystalLength << " Depth " << crystalDepth
0449                               << " crystal name = " << lv->GetName() << " " << dd4hep::dd::noNamespace(lv->GetName())
0450                               << " z of localPoint = " << currentLocalPoint.z() << " take weight = " << weight;
0451 #endif
0452   return weight;
0453 }
0454 
0455 void ECalSD::getBaseNumber(const G4Step* aStep) {
0456   theBaseNumber.reset();
0457   const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0458   int theSize = touch->GetHistoryDepth() + 1;
0459   if (theBaseNumber.getCapacity() < theSize)
0460     theBaseNumber.setSize(theSize);
0461   //Get name and copy numbers
0462   if (theSize > 1) {
0463     for (int ii = 0; ii < theSize; ii++) {
0464       std::string_view name = dd4hep::dd::noNamespace(touch->GetVolume(ii)->GetName());
0465       theBaseNumber.addLevel(std::string(name), touch->GetReplicaNumber(ii));
0466 #ifdef EDM_ML_DEBUG
0467       edm::LogVerbatim("EcalSim") << "ECalSD::getBaseNumber(): Adding level " << ii << ": " << name << "["
0468                                   << touch->GetReplicaNumber(ii) << "]";
0469 #endif
0470     }
0471   }
0472 }
0473 
0474 double ECalSD::getBirkL3(const G4Step* aStep) {
0475   double weight = 1.;
0476   const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0477   double charge = preStepPoint->GetCharge();
0478 
0479   if (charge != 0. && aStep->GetStepLength() > 0.) {
0480     const G4Material* mat = preStepPoint->GetMaterial();
0481     double density = mat->GetDensity();
0482     double dedx = aStep->GetTotalEnergyDeposit() / aStep->GetStepLength();
0483     double rkb = birk1 / density;
0484     if (dedx > 0) {
0485       weight = 1. - birkSlope * log(rkb * dedx);
0486       if (weight < birkCut)
0487         weight = birkCut;
0488       else if (weight > 1.)
0489         weight = 1.;
0490     }
0491 #ifdef EDM_ML_DEBUG
0492     edm::LogVerbatim("EcalSim") << "ECalSD::getBirkL3 in " << dd4hep::dd::noNamespace(mat->GetName()) << " Charge "
0493                                 << charge << " dE/dx " << dedx << " Birk Const " << rkb << " Weight = " << weight
0494                                 << " dE " << aStep->GetTotalEnergyDeposit();
0495 #endif
0496   }
0497   return weight;
0498 }
0499 
0500 bool ECalSD::isXtal(const G4LogicalVolume* lv) { return (xtalLMap.find(lv) != xtalLMap.end()); }
[ 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!