Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-10 02:21:15 

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/Geometry/interface/DD4hep2DDDName.h"
0008 #include "SimG4Core/Notification/interface/TrackInformation.h"
0009 #include "Geometry/EcalCommonData/interface/EcalBarrelNumberingScheme.h"
0010 #include "Geometry/EcalCommonData/interface/EcalBaseNumber.h"
0011 #include "Geometry/EcalCommonData/interface/EcalEndcapNumberingScheme.h"
0012 #include "Geometry/EcalCommonData/interface/EcalPreshowerNumberingScheme.h"
0013 #include "Geometry/EcalCommonData/interface/ESTBNumberingScheme.h"
0014 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0015 #include "DataFormats/Math/interface/GeantUnits.h"
0016 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0017 
0018 #include "FWCore/Framework/interface/ESTransientHandle.h"
0019 #include "FWCore/ServiceRegistry/interface/Service.h"
0020 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0021 
0022 #include "G4LogicalVolumeStore.hh"
0023 #include "G4LogicalVolume.hh"
0024 #include "G4Step.hh"
0025 #include "G4Track.hh"
0026 #include "G4VProcess.hh"
0027 
0028 #include <CLHEP/Units/SystemOfUnits.h>
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   }
0133 #ifdef EDM_ML_DEBUG
0134   edm::LogVerbatim("EcalSim") << "Constructing a ECalSD  with name " << GetName();
0135 #endif
0136   if (useWeight) {
0137     edm::LogVerbatim("EcalSim") << "ECalSD:: Use of Birks law is set to " << useBirk
0138                                 << " with three constants kB = " << birk1 / bunit << ", C1 = " << birk2
0139                                 << ", C2 = " << birk3 << "\n         Use of L3 parametrization " << useBirkL3
0140                                 << " with slope " << birkSlope << " and cut off " << birkCut << "\n"
0141                                 << "         Slope for Light yield is set to " << slopeLY;
0142   } else {
0143     edm::LogVerbatim("EcalSim") << "ECalSD:: energy deposit is not corrected "
0144                                 << " by Birk or light yield curve";
0145   }
0146 
0147   edm::LogVerbatim("EcalSim") << "ECalSD:: Suppression Flag " << suppressHeavy << "\tprotons below "
0148                               << kmaxProton / CLHEP::MeV << " MeV,\tneutrons below " << kmaxNeutron / CLHEP::MeV
0149                               << " MeV,\tions below " << kmaxIon / CLHEP::MeV << " MeV \n\tDepth1 Name = " << depth1Name
0150                               << "\tDepth2 Name = " << depth2Name << "\n\tstoreRL " << storeRL << ":" << scaleRL
0151                               << "\tstoreLayerTimeSim " << storeLayerTimeSim << "\n\ttime Granularity "
0152                               << p.getParameter<edm::ParameterSet>("ECalSD").getParameter<double>("TimeSliceUnit")
0153                               << " ns";
0154   if (useWeight)
0155     initMap();
0156 #ifdef plotDebug
0157   edm::Service<TFileService> tfile;
0158   if (tfile.isAvailable()) {
0159     TFileDirectory ecDir = tfile->mkdir("ProfileFromECalSD");
0160     static const std::string ctype[4] = {"EB", "EBref", "EE", "EERef"};
0161     for (int k = 0; k < 4; ++k) {
0162       std::string name = "ECLL_" + ctype[k];
0163       std::string title = "Local vs Global for " + ctype[k];
0164       double xmin = (k > 1) ? 3000.0 : 1000.0;
0165       g2L_[k] = ecDir.make<TH2F>(name.c_str(), title.c_str(), 100, xmin, xmin + 1000., 100, 0.0, 3000.);
0166     }
0167   } else {
0168     for (int k = 0; k < 4; ++k)
0169       g2L_[k] = 0;
0170   }
0171 #endif
0172   if (dump) {
0173     const auto& lvNames = clg.logicalNames(name);
0174     EcalDumpGeometry geom(lvNames, depth1Name, depth2Name, type);
0175     geom.update();
0176   }
0177 }
0178 
0179 ECalSD::~ECalSD() {}
0180 
0181 double ECalSD::getEnergyDeposit(const G4Step* aStep) {
0182   const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0183   const G4Track* theTrack = aStep->GetTrack();
0184   double edep = aStep->GetTotalEnergyDeposit();
0185 
0186   // take into account light collection curve for crystals
0187   double weight = 1.;
0188   if (suppressHeavy) {
0189     TrackInformation* trkInfo = (TrackInformation*)(theTrack->GetUserInformation());
0190     if (trkInfo) {
0191       int pdg = theTrack->GetDefinition()->GetPDGEncoding();
0192       if (!(trkInfo->isPrimary())) {  // Only secondary particles
0193         double ke = theTrack->GetKineticEnergy();
0194         if (((pdg / 1000000000 == 1 && ((pdg / 10000) % 100) > 0 && ((pdg / 10) % 100) > 0)) && (ke < kmaxIon))
0195           weight = 0;
0196         if ((pdg == 2212) && (ke < kmaxProton))
0197           weight = 0;
0198         if ((pdg == 2112) && (ke < kmaxNeutron))
0199           weight = 0;
0200       }
0201     }
0202   }
0203   const G4LogicalVolume* lv = preStepPoint->GetTouchable()->GetVolume(0)->GetLogicalVolume();
0204   double wt1 = 1.0;
0205   if (useWeight && !any(noWeight, lv)) {
0206     weight *= curve_LY(lv);
0207     if (useBirk) {
0208       if (useBirkL3)
0209         weight *= getBirkL3(aStep);
0210       else
0211         weight *= getAttenuation(aStep, birk1, birk2, birk3);
0212     }
0213     wt1 = getResponseWt(theTrack);
0214   }
0215   edep *= weight * wt1;
0216   // Russian Roulette
0217   double wt2 = theTrack->GetWeight();
0218   if (wt2 > 0.0) {
0219     edep *= wt2;
0220   }
0221 #ifdef EDM_ML_DEBUG
0222   edm::LogVerbatim("EcalSim") << lv->GetName() << " "
0223                               << DD4hep2DDDName::noNameSpace(static_cast<std::string>(lv->GetName()))
0224                               << " Light Collection Efficiency " << weight << ":" << wt1 << " wt2= " << wt2
0225                               << " Weighted Energy Deposit " << edep / CLHEP::MeV << " MeV at "
0226                               << preStepPoint->GetPosition();
0227 #endif
0228   return edep;
0229 }
0230 
0231 double ECalSD::EnergyCorrected(const G4Step& step, const G4Track* track) {
0232   const G4StepPoint* hitPoint = step.GetPreStepPoint();
0233   const G4LogicalVolume* lv = hitPoint->GetTouchable()->GetVolume(0)->GetLogicalVolume();
0234   if (lv->GetSensitiveDetector() != this)
0235     return 0.0;
0236 
0237   double edep = step.GetTotalEnergyDeposit();
0238   if (useWeight && !any(noWeight, lv)) {
0239     currentLocalPoint = setToLocal(hitPoint->GetPosition(), hitPoint->GetTouchable());
0240     auto ite = xtalLMap.find(lv);
0241     crystalLength = (ite == xtalLMap.end()) ? 230._mm : std::abs(ite->second);
0242     crystalDepth = (ite == xtalLMap.end()) ? 0.0 : (std::abs(0.5 * (ite->second) + currentLocalPoint.z()));
0243     edep *= curve_LY(lv) * getResponseWt(track);
0244   }
0245   return edep;
0246 }
0247 
0248 int ECalSD::getTrackID(const G4Track* aTrack) {
0249   int primaryID(0);
0250   if (storeTrack && depth > 0) {
0251     forceSave = true;
0252     primaryID = aTrack->GetTrackID();
0253   } else {
0254     primaryID = CaloSD::getTrackID(aTrack);
0255   }
0256   return primaryID;
0257 }
0258 
0259 uint16_t ECalSD::getDepth(const G4Step* aStep) {
0260   // this method should be called first at a step
0261   const G4StepPoint* hitPoint = aStep->GetPreStepPoint();
0262   currentLocalPoint = setToLocal(hitPoint->GetPosition(), hitPoint->GetTouchable());
0263   const G4LogicalVolume* lv = hitPoint->GetTouchable()->GetVolume(0)->GetLogicalVolume();
0264 
0265   auto ite = xtalLMap.find(lv);
0266   crystalLength = (ite == xtalLMap.end()) ? 230._mm : std::abs(ite->second);
0267   crystalDepth = (ite == xtalLMap.end()) ? 0.0 : (std::abs(0.5 * (ite->second) + currentLocalPoint.z()));
0268   depth = any(useDepth1, lv) ? 1 : (any(useDepth2, lv) ? 2 : 0);
0269   uint16_t depth1(0), depth2(0);
0270   if (storeRL) {
0271     depth1 = (ite == xtalLMap.end()) ? 0 : (((ite->second) >= 0) ? 0 : PCaloHit::kEcalDepthRefz);
0272     depth2 = getRadiationLength(hitPoint, lv);
0273     depth |= (((depth2 & PCaloHit::kEcalDepthMask) << PCaloHit::kEcalDepthOffset) | depth1);
0274   } else if (storeLayerTimeSim) {
0275     depth2 = getLayerIDForTimeSim();
0276     depth |= ((depth2 & PCaloHit::kEcalDepthMask) << PCaloHit::kEcalDepthOffset);
0277   }
0278 #ifdef EDM_ML_DEBUG
0279   if (isXtal(lv))
0280     edm::LogVerbatim("EcalSimX") << "ECalSD::Volume " << lv->GetName() << " DetId " << std::hex << setDetUnitId(aStep)
0281                                  << std::dec << " Global " << (hitPoint->GetPosition()).rho() << ":"
0282                                  << (hitPoint->GetPosition()).z() << " Local Z " << currentLocalPoint.z() << " Depth "
0283                                  << crystalDepth;
0284   edm::LogVerbatim("EcalSim") << "ECalSD::Depth " << std::hex << depth1 << ":" << depth2 << ":" << depth << std::dec
0285                               << " L " << (ite == xtalLMap.end()) << ":" << ite->second << " local "
0286                               << currentLocalPoint << " Crystal length " << crystalLength << ":" << crystalDepth;
0287 #endif
0288   return depth;
0289 }
0290 
0291 uint16_t ECalSD::getRadiationLength(const G4StepPoint* hitPoint, const G4LogicalVolume* lv) {
0292   uint16_t thisX0 = 0;
0293   if (useWeight) {
0294     double radl = hitPoint->GetMaterial()->GetRadlen();
0295     thisX0 = (uint16_t)floor(scaleRL * crystalDepth / radl);
0296 #ifdef plotDebug
0297     const std::string lvname = DD4hep2DDDName::noNameSpace(static_cast<std::string>(lv->GetName()));
0298     int k1 = (lvname.find("EFRY") != std::string::npos) ? 2 : 0;
0299     int k2 = (lvname.find("refl") != std::string::npos) ? 1 : 0;
0300     int kk = k1 + k2;
0301     double rz = (k1 == 0) ? (hitPoint->GetPosition()).rho() : std::abs((hitPoint->GetPosition()).z());
0302     edm::LogVerbatim("EcalSim") << lvname << " # " << k1 << ":" << k2 << ":" << kk << " rz " << rz << " D " << thisX0;
0303     g2L_[kk]->Fill(rz, thisX0);
0304 #endif
0305 #ifdef EDM_ML_DEBUG
0306     G4ThreeVector localPoint = setToLocal(hitPoint->GetPosition(), hitPoint->GetTouchable());
0307     edm::LogVerbatim("EcalSim") << lv->GetName() << " "
0308                                 << DD4hep2DDDName::noNameSpace(static_cast<std::string>(lv->GetName())) << " Global "
0309                                 << hitPoint->GetPosition() << ":" << (hitPoint->GetPosition()).rho() << " Local "
0310                                 << localPoint << " Crystal Length " << crystalLength << " Radl " << radl
0311                                 << " crystalDepth " << crystalDepth << " Index " << thisX0 << " : "
0312                                 << getLayerIDForTimeSim();
0313 #endif
0314   }
0315   return thisX0;
0316 }
0317 
0318 uint16_t ECalSD::getLayerIDForTimeSim() {
0319   const double invLayerSize = 0.1;  //layer size in 1/mm
0320   return (int)crystalDepth * invLayerSize;
0321 }
0322 
0323 uint32_t ECalSD::setDetUnitId(const G4Step* aStep) {
0324   if (numberingScheme_.get() == nullptr) {
0325     return EBDetId(1, 1)();
0326   } else {
0327     getBaseNumber(aStep);
0328     return numberingScheme_->getUnitID(theBaseNumber);
0329   }
0330 }
0331 
0332 void ECalSD::setNumberingScheme(EcalNumberingScheme* scheme) {
0333   if (scheme != nullptr) {
0334     edm::LogVerbatim("EcalSim") << "EcalSD: updates numbering scheme for " << GetName();
0335     numberingScheme_.reset(scheme);
0336   }
0337 }
0338 
0339 void ECalSD::initMap() {
0340   std::vector<const G4LogicalVolume*> lvused;
0341   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0342   std::map<const std::string, const G4LogicalVolume*> nameMap;
0343   for (auto lvi = lvs->begin(), lve = lvs->end(); lvi != lve; ++lvi)
0344     nameMap.emplace(DD4hep2DDDName::noNameSpace(static_cast<std::string>((*lvi)->GetName())), *lvi);
0345 
0346   for (unsigned int it = 0; it < ecalSimParameters_->lvNames_.size(); ++it) {
0347     const std::string& matname = ecalSimParameters_->matNames_[it];
0348     const std::string& lvname = ecalSimParameters_->lvNames_[it];
0349     const G4LogicalVolume* lv = nameMap[lvname];
0350     int ibec = (lvname.find("EFRY") == std::string::npos) ? 0 : 1;
0351     int iref = (lvname.find("refl") == std::string::npos) ? 0 : 1;
0352     int type = (ibec + iref == 1) ? 1 : -1;
0353     if (depth1Name != " ") {
0354       if (strncmp(lvname.c_str(), depth1Name.c_str(), 4) == 0) {
0355         if (!any(useDepth1, lv)) {
0356           useDepth1.push_back(lv);
0357 #ifdef EDM_ML_DEBUG
0358           edm::LogVerbatim("EcalSim") << "ECalSD::initMap Logical Volume " << lvname << " in Depth 1 volume list";
0359 #endif
0360         }
0361         const G4LogicalVolume* lvr = nameMap[lvname + "_refl"];
0362         if (lvr != nullptr && !any(useDepth1, lvr)) {
0363           useDepth1.push_back(lvr);
0364 #ifdef EDM_ML_DEBUG
0365           edm::LogVerbatim("EcalSim") << "ECalSD::initMap Logical Volume " << lvname << "_refl"
0366                                       << " in Depth 1 volume list";
0367 #endif
0368         }
0369       }
0370     }
0371     if (depth2Name != " ") {
0372       if (strncmp(lvname.c_str(), depth2Name.c_str(), 4) == 0) {
0373         if (!any(useDepth2, lv)) {
0374           useDepth2.push_back(lv);
0375 #ifdef EDM_ML_DEBUG
0376           edm::LogVerbatim("EcalSim") << "ECalSD::initMap Logical Volume " << lvname << " in Depth 2 volume list";
0377 #endif
0378         }
0379         const G4LogicalVolume* lvr = nameMap[lvname + "_refl"];
0380         if (lvr != nullptr && !any(useDepth2, lvr)) {
0381           useDepth2.push_back(lvr);
0382 #ifdef EDM_ML_DEBUG
0383           edm::LogVerbatim("EcalSim") << "ECalSD::initMap Logical Volume " << lvname << "_refl"
0384                                       << " in Depth 2 volume list";
0385 #endif
0386         }
0387       }
0388     }
0389     if (lv != nullptr) {
0390       if (crystalMat.size() == matname.size() && !strcmp(crystalMat.c_str(), matname.c_str())) {
0391         if (!any(lvused, lv)) {
0392           lvused.push_back(lv);
0393           double dz = ecalSimParameters_->dzs_[it];
0394           xtalLMap.insert(std::pair<const G4LogicalVolume*, double>(lv, dz * type));
0395           lv = nameMap[lvname + "_refl"];
0396           if (lv != nullptr) {
0397             xtalLMap.insert(std::pair<const G4LogicalVolume*, double>(lv, -dz * type));
0398           }
0399         }
0400       } else {
0401         if (!any(noWeight, lv)) {
0402           noWeight.push_back(lv);
0403 #ifdef EDM_ML_DEBUG
0404           edm::LogVerbatim("EcalSim") << "ECalSD::initMap Logical Volume " << lvname << " Material " << matname
0405                                       << " in noWeight list";
0406 #endif
0407         }
0408         lv = nameMap[lvname];
0409         if (lv != nullptr && !any(noWeight, lv)) {
0410           noWeight.push_back(lv);
0411 #ifdef EDM_ML_DEBUG
0412           edm::LogVerbatim("EcalSim") << "ECalSD::initMap Logical Volume " << lvname << " Material " << matname
0413                                       << " in noWeight list";
0414 #endif
0415         }
0416       }
0417     }
0418   }
0419 #ifdef EDM_ML_DEBUG
0420   edm::LogVerbatim("EcalSim") << "ECalSD: Length Table:";
0421   int i = 0;
0422   for (auto ite : xtalLMap) {
0423     std::string name("Unknown");
0424     if (ite.first != nullptr)
0425       name = DD4hep2DDDName::noNameSpace(static_cast<std::string>((ite.first)->GetName()));
0426     edm::LogVerbatim("EcalSim") << " " << i << " " << ite.first << " " << name << " L = " << ite.second;
0427     ++i;
0428   }
0429 #endif
0430 }
0431 
0432 double ECalSD::curve_LY(const G4LogicalVolume* lv) {
0433   double weight = 1.;
0434   if (ageingWithSlopeLY) {
0435     //position along the crystal in mm from 0 to 230 (in EB)
0436     if (crystalDepth >= -0.1 || crystalDepth <= crystalLength + 0.1)
0437       weight = ageing.calcLightCollectionEfficiencyWeighted(currentID[0].unitID(), crystalDepth / crystalLength);
0438   } else {
0439     double dapd = crystalLength - crystalDepth;
0440     if (dapd >= -0.1 || dapd <= crystalLength + 0.1) {
0441       if (dapd <= 100.)
0442         weight = 1.0 + slopeLY - dapd * 0.01 * slopeLY;
0443     } else {
0444       edm::LogWarning("EcalSim") << "ECalSD: light coll curve : wrong distance "
0445                                  << "to APD " << dapd << " crlength = " << crystalLength << ":" << crystalDepth
0446                                  << " crystal name = " << lv->GetName() << " "
0447                                  << DD4hep2DDDName::noNameSpace(static_cast<std::string>(lv->GetName()))
0448                                  << " z of localPoint = " << currentLocalPoint.z() << " take weight = " << weight;
0449     }
0450   }
0451 #ifdef EDM_ML_DEBUG
0452   edm::LogVerbatim("EcalSim") << "ECalSD: light coll curve : crlength = " << crystalLength << " Depth " << crystalDepth
0453                               << " crystal name = " << lv->GetName() << " "
0454                               << DD4hep2DDDName::noNameSpace(static_cast<std::string>(lv->GetName()))
0455                               << " z of localPoint = " << currentLocalPoint.z() << " take weight = " << weight;
0456 #endif
0457   return weight;
0458 }
0459 
0460 void ECalSD::getBaseNumber(const G4Step* aStep) {
0461   theBaseNumber.reset();
0462   const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0463   int theSize = touch->GetHistoryDepth() + 1;
0464   if (theBaseNumber.getCapacity() < theSize)
0465     theBaseNumber.setSize(theSize);
0466   //Get name and copy numbers
0467   if (theSize > 1) {
0468     for (int ii = 0; ii < theSize; ii++) {
0469       std::string name = DD4hep2DDDName::noNameSpace(static_cast<std::string>(touch->GetVolume(ii)->GetName()));
0470       theBaseNumber.addLevel(name, touch->GetReplicaNumber(ii));
0471 #ifdef EDM_ML_DEBUG
0472       edm::LogVerbatim("EcalSim") << "ECalSD::getBaseNumber(): Adding level " << ii << ": " << name << "["
0473                                   << touch->GetReplicaNumber(ii) << "]";
0474 #endif
0475     }
0476   }
0477 }
0478 
0479 double ECalSD::getBirkL3(const G4Step* aStep) {
0480   double weight = 1.;
0481   const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0482   double charge = preStepPoint->GetCharge();
0483 
0484   if (charge != 0. && aStep->GetStepLength() > 0.) {
0485     const G4Material* mat = preStepPoint->GetMaterial();
0486     double density = mat->GetDensity();
0487     double dedx = aStep->GetTotalEnergyDeposit() / aStep->GetStepLength();
0488     double rkb = birk1 / density;
0489     if (dedx > 0) {
0490       weight = 1. - birkSlope * log(rkb * dedx);
0491       if (weight < birkCut)
0492         weight = birkCut;
0493       else if (weight > 1.)
0494         weight = 1.;
0495     }
0496 #ifdef EDM_ML_DEBUG
0497     edm::LogVerbatim("EcalSim") << "ECalSD::getBirkL3 in "
0498                                 << DD4hep2DDDName::noNameSpace(static_cast<std::string>(mat->GetName())) << " Charge "
0499                                 << charge << " dE/dx " << dedx << " Birk Const " << rkb << " Weight = " << weight
0500                                 << " dE " << aStep->GetTotalEnergyDeposit();
0501 #endif
0502   }
0503   return weight;
0504 }
0505 
0506 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!