Forward/src/MtdSD.cc

0001 //#define EDM_ML_DEBUG
0002
0003 #include "SimG4CMS/Forward/interface/MtdSD.h"
0004
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007
0008 #include "Geometry/MTDCommonData/interface/BTLNumberingScheme.h"
0009 #include "Geometry/MTDCommonData/interface/ETLNumberingScheme.h"
0010 #include "DataFormats/ForwardDetId/interface/MTDDetId.h"
0011 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0012
0013 #include "G4Track.hh"
0014 #include "G4Step.hh"
0015 #include "G4StepPoint.hh"
0016
0017 #include <iostream>
0018
0019 using namespace MtdHitCategory;
0020
0021 //-------------------------------------------------------------------
0022 MtdSD::MtdSD(const std::string& name,
0023              const SensitiveDetectorCatalog& clg,
0024              edm::ParameterSet const& p,
0025              const SimTrackManager* manager)
0026     : TimingSD(name, clg, manager), numberingScheme(nullptr) {
0027   //Parameters
0028   edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("MtdSD");
0029   int verbn = m_p.getUntrackedParameter<int>("Verbosity");
0030
0031   SetVerboseLevel(verbn);
0032
0033   MTDNumberingScheme* scheme = nullptr;
0034   if (name == "FastTimerHitsBarrel") {
0035     scheme = dynamic_cast<MTDNumberingScheme*>(new BTLNumberingScheme());
0036     isBTL = true;
0037   } else if (name == "FastTimerHitsEndcap") {
0038     scheme = dynamic_cast<MTDNumberingScheme*>(new ETLNumberingScheme());
0039     isETL = true;
0040   } else {
0041     scheme = nullptr;
0042     edm::LogWarning("MtdSim") << "MtdSD: ReadoutName not supported";
0043   }
0044   if (scheme)
0045     setNumberingScheme(scheme);
0046
0047   energyCut = m_p.getParameter<double>("EnergyThresholdForPersistencyInGeV") * CLHEP::GeV;     //default must be 0.5
0048   energyHistoryCut = m_p.getParameter<double>("EnergyThresholdForHistoryInGeV") * CLHEP::GeV;  //default must be 0.05
0049
0050   setCuts(energyCut, energyHistoryCut);
0051
0052   double newTimeFactor = 1. / m_p.getParameter<double>("TimeSliceUnit");
0053   edm::LogVerbatim("MtdSim") << "New time factor = " << newTimeFactor;
0054   setTimeFactor(newTimeFactor);
0055
0056   edm::LogVerbatim("MtdSim") << "MtdSD: Instantiation completed for " << name;
0057 }
0058
0059 MtdSD::~MtdSD() {}
0060
0061 uint32_t MtdSD::setDetUnitId(const G4Step* aStep) {
0062   if (numberingScheme == nullptr) {
0063     return MTDDetId();
0064   } else {
0065     getBaseNumber(aStep);
0066 #ifdef EDM_ML_DEBUG
0067     edm::LogVerbatim("MtdSim") << "DetId = " << numberingScheme->getUnitID(theBaseNumber);
0068 #endif
0069     return numberingScheme->getUnitID(theBaseNumber);
0070   }
0071 }
0072
0073 void MtdSD::setNumberingScheme(MTDNumberingScheme* scheme) {
0074   if (scheme != nullptr) {
0075     edm::LogVerbatim("MtdSim") << "MtdSD: updates numbering scheme for " << GetName();
0076     if (numberingScheme)
0077       delete numberingScheme;
0078     numberingScheme = scheme;
0079   }
0080 }
0081
0082 void MtdSD::getBaseNumber(const G4Step* aStep) {
0083   theBaseNumber.reset();
0084   const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0085   int theSize = touch->GetHistoryDepth() + 1;
0086   if (theBaseNumber.getCapacity() < theSize)
0087     theBaseNumber.setSize(theSize);
0088   //Get name and copy numbers
0089   if (theSize > 1) {
0090 #ifdef EDM_ML_DEBUG
0091     edm::LogVerbatim("MtdSim") << "Building MTD basenumber:";
0092 #endif
0093     for (int ii = 0; ii < theSize; ii++) {
0094       theBaseNumber.addLevel(touch->GetVolume(ii)->GetName(), touch->GetReplicaNumber(ii));
0095 #ifdef EDM_ML_DEBUG
0096       edm::LogVerbatim("MtdSim") << "MtdSD::getBaseNumber(): Adding level " << ii << ": "
0097                                  << touch->GetVolume(ii)->GetName() << "[" << touch->GetReplicaNumber(ii) << "]";
0098 #endif
0099     }
0100   }
0101 }
0102
0103 int MtdSD::getTrackID(const G4Track* aTrack) {
0104   int theID = aTrack->GetTrackID();
0105   TrackInformation* trkInfo = cmsTrackInformation(aTrack);
0106   const G4String& rname = aTrack->GetVolume()->GetLogicalVolume()->GetRegion()->GetName();
0107   if (trkInfo != nullptr) {
0108 #ifdef EDM_ML_DEBUG
0109     trkInfo->Print();
0110 #endif
0111     if (!trkInfo->storeTrack()) {
0112       theID = trkInfo->idLastStoredAncestor();
0113     }
0114     if (rname == "FastTimerRegionSensBTL") {
0115       if (trkInfo->isInTrkFromBackscattering()) {
0116         theID = PSimHit::addTrackIdOffset(theID, k_idFromCaloOffset);
0117       } else if (trkInfo->isExtSecondary() && !trkInfo->isInTrkFromBackscattering() && !trkInfo->storeTrack()) {
0118         theID = PSimHit::addTrackIdOffset(theID, k_idsecOffset);
0119       } else if (trkInfo->isBTLlooper()) {
0120         theID = PSimHit::addTrackIdOffset(theID, k_idloopOffset);
0121       }
0122 #ifdef EDM_ML_DEBUG
0123       edm::LogVerbatim("MtdSim") << "MtdSD: Track ID: " << aTrack->GetTrackID()
0124                                  << " BTL Track ID: " << trkInfo->mcTruthID() << ":" << theID;
0125 #endif
0126     } else if (rname == "FastTimerRegionSensETL") {
0127       if (hitClassID == k_idETLfromBack) {
0128         theID = PSimHit::addTrackIdOffset(theID, k_idETLfromBack);
0129       }
0130 #ifdef EDM_ML_DEBUG
0131       edm::LogVerbatim("MtdSim") << "MtdSD: Track ID: " << aTrack->GetTrackID()
0132                                  << " ETL Track ID: " << trkInfo->mcTruthID() << ":" << theID;
0133 #endif
0134       // In the case of ECAL GFlash fast spot may be inside MTD and should be ignored
0135     } else if (rname == "EcalRegion") {
0136       theID = -2;
0137     } else {
0138       throw cms::Exception("MtdSDError") << "MtdSD called in incorrect region " << rname;
0139     }
0140   } else {
0141 #ifdef EDM_ML_DEBUG
0142     edm::LogWarning("MtdSim") << "MtdSD: Problem with primaryID **** set by force to TkID **** " << theID;
0143 #endif
0144   }
0145   return theID;
0146 }
0147
0148 void MtdSD::setHitClassID(const G4Step* aStep) {
0149   TrackInformation* trkInfo = cmsTrackInformation(aStep->GetTrack());
0150   const G4String& rname = aStep->GetTrack()->GetVolume()->GetLogicalVolume()->GetRegion()->GetName();
0151   if (rname == "FastTimerRegionSensETL") {
0152     double zin = std::abs(aStep->GetPreStepPoint()->GetPosition().z());
0153     double zout = std::abs(aStep->GetPostStepPoint()->GetPosition().z());
0154     if (zout - zin < 0.) {
0155       hitClassID = k_idETLfromBack;
0156       trkInfo->setETLfromBack();
0157     } else {
0158       hitClassID = 0;
0159       trkInfo->setETLfromFront();
0160     }
0161   }
0162 }