Muon/src/MuonSensitiveDetector.cc

0001 #include "SimG4CMS/Muon/interface/MuonSensitiveDetector.h"
0002 #include "SimG4CMS/Muon/interface/MuonSlaveSD.h"
0003 #include "SimG4CMS/Muon/interface/MuonEndcapFrameRotation.h"
0004 #include "SimG4CMS/Muon/interface/MuonRPCFrameRotation.h"
0005 #include "SimG4CMS/Muon/interface/MuonGEMFrameRotation.h"
0006 #include "SimG4CMS/Muon/interface/MuonME0FrameRotation.h"
0007 #include "Geometry/MuonNumbering/interface/MuonSubDetector.h"
0008
0009 #include "SimG4CMS/Muon/interface/SimHitPrinter.h"
0010 #include "SimDataFormats/TrackingHit/interface/UpdatablePSimHit.h"
0011 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0012 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0013 #include "DataFormats/MuonDetId/interface/GEMDetId.h"
0014 #include "DataFormats/MuonDetId/interface/ME0DetId.h"
0015 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0016
0017 #include "SimG4CMS/Muon/interface/MuonG4Numbering.h"
0018 #include "Geometry/MuonNumbering/interface/MuonGeometryConstants.h"
0019 #include "Geometry/MuonNumbering/interface/MuonBaseNumber.h"
0020 #include "Geometry/MuonNumbering/interface/MuonSimHitNumberingScheme.h"
0021
0022 #include "SimG4Core/Notification/interface/TrackInformation.h"
0023 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0024 #include "SimG4Core/Physics/interface/G4ProcessTypeEnumerator.h"
0025
0026 #include "G4VProcess.hh"
0027 #include "G4EventManager.hh"
0028 #include <CLHEP/Units/SystemOfUnits.h>
0029 #include "G4Step.hh"
0030 #include "G4StepPoint.hh"
0031 #include "G4Track.hh"
0032
0033 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0034
0035 #include <iostream>
0036 #include <sstream>
0037
0038 //#define EDM_ML_DEBUG
0039
0040 MuonSensitiveDetector::MuonSensitiveDetector(const std::string& name,
0041                                              const MuonOffsetMap* offmap,
0042                                              const MuonGeometryConstants& constants,
0043                                              const SensitiveDetectorCatalog& clg,
0044                                              edm::ParameterSet const& p,
0045                                              const SimTrackManager* manager)
0046     : SensitiveTkDetector(name, clg),
0047       thePV(nullptr),
0048       theHit(nullptr),
0049       theDetUnitId(0),
0050       newDetUnitId(0),
0051       theTrackID(0),
0052       thePrinter(nullptr),
0053       theManager(manager) {
0054   // Here simply create 1 MuonSlaveSD for the moment
0055   //
0056   bool dd4hep = p.getParameter<bool>("g4GeometryDD4hepSource");
0057   edm::ParameterSet muonSD = p.getParameter<edm::ParameterSet>("MuonSD");
0058   printHits_ = muonSD.getParameter<bool>("PrintHits");
0059   ePersistentCutGeV_ = muonSD.getParameter<double>("EnergyThresholdForPersistency") / CLHEP::GeV;  //Default 1. GeV
0060   allMuonsPersistent_ = muonSD.getParameter<bool>("AllMuonsPersistent");
0061   haveDemo_ = muonSD.getParameter<bool>("HaveDemoChambers");
0062   demoGEM_ = muonSD.getParameter<bool>("UseDemoHitGEM");
0063   demoRPC_ = muonSD.getParameter<bool>("UseDemoHitRPC");
0064   removeGEMHits_ = muonSD.getParameter<std::vector<int32_t>>("RemoveGEMHits");
0065
0066 #ifdef EDM_ML_DEBUG
0067   edm::LogVerbatim("MuonSim") << "create MuonSubDetector " << name << " with dd4hep flag " << dd4hep
0068                               << " Flags for Demonstration chambers " << haveDemo_ << " for GEM " << demoGEM_
0069                               << " for RPC " << demoRPC_ << " " << removeGEMHits_.size() << " GEM Hits to be removed";
0070   if (!removeGEMHits_.empty()) {
0071     std::ostringstream st1;
0072     for (const auto& id : removeGEMHits_)
0073       st1 << " " << id;
0074     edm::LogVerbatim("MuonSim") << "IDs to be removed " << st1.str();
0075   }
0076 #endif
0077   detector = new MuonSubDetector(name);
0078
0079   G4String sdet = "unknown";
0080   if (detector->isEndcap()) {
0081     theRotation = new MuonEndcapFrameRotation();
0082     sdet = "Endcap";
0083   } else if (detector->isRPC()) {
0084     theRotation = new MuonRPCFrameRotation(constants, offmap, dd4hep);
0085     sdet = "RPC";
0086   } else if (detector->isGEM()) {
0087     theRotation = new MuonGEMFrameRotation(constants);
0088     sdet = "GEM";
0089   } else if (detector->isME0()) {
0090     theRotation = new MuonME0FrameRotation(constants);
0091     sdet = "ME0";
0092   } else {
0093     theRotation = new MuonFrameRotation();
0094     if (detector->isBarrel())
0095       sdet = "Barrel";
0096   }
0097   slaveMuon = new MuonSlaveSD(detector, theManager);
0098   numbering = new MuonSimHitNumberingScheme(detector, constants);
0099   g4numbering = new MuonG4Numbering(constants, offmap, dd4hep);
0100
0101   if (printHits_) {
0102     thePrinter = new SimHitPrinter("HitPositionOSCAR.dat");
0103   }
0104
0105   edm::LogVerbatim("MuonSim") << " of type " << sdet << " <" << GetName() << "> EnergyThresholdForPersistency(GeV) "
0106                               << ePersistentCutGeV_ / CLHEP::GeV << " allMuonsPersistent: " << allMuonsPersistent_;
0107
0108   theG4ProcessTypeEnumerator = new G4ProcessTypeEnumerator;
0109 }
0110
0111 MuonSensitiveDetector::~MuonSensitiveDetector() {
0112   delete g4numbering;
0113   delete numbering;
0114   delete slaveMuon;
0115   delete theRotation;
0116   delete detector;
0117   delete theG4ProcessTypeEnumerator;
0118 }
0119
0120 void MuonSensitiveDetector::update(const BeginOfEvent* i) {
0121   clearHits();
0122   //----- Initialize variables to check if two steps belong to same hit
0123   thePV = nullptr;
0124   theDetUnitId = 0;
0125   theTrackID = 0;
0126 }
0127
0128 void MuonSensitiveDetector::clearHits() {
0129 #ifdef EDM_ML_DEBUG
0130   edm::LogVerbatim("MuonSim") << "MuonSensitiveDetector::clearHits";
0131 #endif
0132   slaveMuon->Initialize();
0133 }
0134
0135 bool MuonSensitiveDetector::ProcessHits(G4Step* aStep, G4TouchableHistory* ROhist) {
0136 #ifdef EDM_ML_DEBUG
0137   edm::LogVerbatim("MuonSim") << " MuonSensitiveDetector::ProcessHits " << InitialStepPosition(aStep, WorldCoordinates);
0138 #endif
0139
0140   if (aStep->GetTotalEnergyDeposit() > 0.) {
0141     newDetUnitId = setDetUnitId(aStep);
0142 #ifdef EDM_ML_DEBUG
0143     G4VPhysicalVolume* vol = aStep->GetPreStepPoint()->GetTouchable()->GetVolume(0);
0144     std::string namx = static_cast<std::string>(vol->GetName());
0145     std::string name = namx.substr(0, 2);
0146     if (name == "RE")
0147       edm::LogVerbatim("MuonSim") << "DETID " << namx << " " << RPCDetId(newDetUnitId);
0148 #endif
0149     if (newHit(aStep)) {
0150       saveHit();
0151       createHit(aStep);
0152     } else {
0153       updateHit(aStep);
0154     }
0155     thePV = aStep->GetPreStepPoint()->GetPhysicalVolume();
0156     theTrackID = aStep->GetTrack()->GetTrackID();
0157     theDetUnitId = newDetUnitId;
0158   }
0159   return true;
0160 }
0161
0162 uint32_t MuonSensitiveDetector::setDetUnitId(const G4Step* aStep) {
0163   MuonBaseNumber num = g4numbering->PhysicalVolumeToBaseNumber(aStep);
0164
0165 #ifdef EDM_ML_DEBUG
0166   std::stringstream MuonBaseNumber;
0167   MuonBaseNumber << "MuonNumbering :: number of levels = " << num.getLevels() << std::endl;
0168   MuonBaseNumber << "Level \t SuperNo \t BaseNo" << std::endl;
0169   for (int level = 1; level <= num.getLevels(); level++) {
0170     MuonBaseNumber << level << " \t " << num.getSuperNo(level) << " \t " << num.getBaseNo(level) << std::endl;
0171   }
0172   std::string MuonBaseNumbr = MuonBaseNumber.str();
0173
0174   edm::LogVerbatim("MuonSim") << "MuonSensitiveDetector::setDetUnitId :: " << MuonBaseNumbr;
0175   edm::LogVerbatim("MuonSim") << "MuonSensitiveDetector::setDetUnitId :: MuonDetUnitId = "
0176                               << (numbering->baseNumberToUnitNumber(num));
0177 #endif
0178   return numbering->baseNumberToUnitNumber(num);
0179 }
0180
0181 bool MuonSensitiveDetector::newHit(const G4Step* aStep) {
0182   return (!theHit || (aStep->GetTrack()->GetTrackID() != theTrackID) ||
0183           (aStep->GetPreStepPoint()->GetPhysicalVolume() != thePV) || newDetUnitId != theDetUnitId);
0184 }
0185
0186 void MuonSensitiveDetector::createHit(const G4Step* aStep) {
0187   Local3DPoint theEntryPoint;
0188   Local3DPoint theExitPoint;
0189
0190   if (detector->isBarrel()) {
0191     // 1 levels up
0192     theEntryPoint = cmsUnits(theRotation->transformPoint(InitialStepPositionVsParent(aStep, 1), aStep));
0193     theExitPoint = cmsUnits(theRotation->transformPoint(FinalStepPositionVsParent(aStep, 1), aStep));
0194   } else if (detector->isEndcap()) {
0195     // save local z at current level
0196     theEntryPoint = theRotation->transformPoint(InitialStepPosition(aStep, LocalCoordinates), aStep);
0197     theExitPoint = theRotation->transformPoint(FinalStepPosition(aStep, LocalCoordinates), aStep);
0198     float zentry = theEntryPoint.z();
0199     float zexit = theExitPoint.z();
0200     // 4 levels up
0201     Local3DPoint tempEntry = theRotation->transformPoint(InitialStepPositionVsParent(aStep, 4), aStep);
0202     Local3DPoint tempExit = theRotation->transformPoint(FinalStepPositionVsParent(aStep, 4), aStep);
0203     // reset local z from z wrt deep-parent volume to z wrt low-level volume
0204     theEntryPoint = cmsUnits(Local3DPoint(tempEntry.x(), tempEntry.y(), zentry));
0205     theExitPoint = cmsUnits(Local3DPoint(tempExit.x(), tempExit.y(), zexit));
0206   } else {
0207     theEntryPoint = cmsUnits(theRotation->transformPoint(InitialStepPosition(aStep, LocalCoordinates), aStep));
0208     theExitPoint = cmsUnits(theRotation->transformPoint(FinalStepPosition(aStep, LocalCoordinates), aStep));
0209   }
0210
0211   const G4Track* theTrack = aStep->GetTrack();
0212   const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0213
0214   float thePabs = preStepPoint->GetMomentum().mag() / CLHEP::GeV;
0215   float theTof = preStepPoint->GetGlobalTime() / CLHEP::nanosecond;
0216   float theEnergyLoss = aStep->GetTotalEnergyDeposit() / CLHEP::GeV;
0217   int theParticleType = G4TrackToParticleID::particleID(theTrack);
0218
0219   theDetUnitId = newDetUnitId;
0220   thePV = preStepPoint->GetPhysicalVolume();
0221   theTrackID = theTrack->GetTrackID();
0222
0223   // convert momentum direction it to local frame
0224   const G4ThreeVector& gmd = preStepPoint->GetMomentumDirection();
0225   G4ThreeVector lmd = static_cast<const G4TouchableHistory*>(preStepPoint->GetTouchable())
0226                           ->GetHistory()
0227                           ->GetTopTransform()
0228                           .TransformAxis(gmd);
0229   Local3DPoint lnmd = ConvertToLocal3DPoint(lmd);
0230   lnmd = theRotation->transformPoint(lnmd, aStep);
0231   float theThetaAtEntry = lnmd.theta();
0232   float thePhiAtEntry = lnmd.phi();
0233
0234   theHit = new UpdatablePSimHit(theEntryPoint,
0235                                 theExitPoint,
0236                                 thePabs,
0237                                 theTof,
0238                                 theEnergyLoss,
0239                                 theParticleType,
0240                                 theDetUnitId,
0241                                 theTrackID,
0242                                 theThetaAtEntry,
0243                                 thePhiAtEntry,
0244                                 theG4ProcessTypeEnumerator->processId(theTrack->GetCreatorProcess()));
0245
0246   // Make track persistent
0247   int thePID = std::abs(theTrack->GetDefinition()->GetPDGEncoding());
0248   //---VI - in parameters cut in energy is declared but applied to momentum
0249   if (thePabs > ePersistentCutGeV_ || (thePID == 13 && allMuonsPersistent_)) {
0250     TrackInformation* info = cmsTrackInformation(theTrack);
0251     info->setStoreTrack();
0252     info->setIdLastStoredAncestor(theTrack->GetTrackID());
0253   }
0254
0255 #ifdef EDM_ML_DEBUG
0256   edm::LogVerbatim("MuonSim") << "=== NEW Muon hit for " << GetName() << " Edep(GeV)= " << theEnergyLoss << " "
0257                               << thePV->GetLogicalVolume()->GetName();
0258   const G4VProcess* p = aStep->GetPostStepPoint()->GetProcessDefinedStep();
0259   const G4VProcess* p2 = aStep->GetPreStepPoint()->GetProcessDefinedStep();
0260   G4String sss = "";
0261   if (p)
0262     sss += " POST PROCESS: " + p->GetProcessName();
0263   if (p2)
0264     sss += ";  PRE  PROCESS: " + p2->GetProcessName();
0265   if (!sss.empty())
0266     edm::LogVerbatim("MuonSim") << sss;
0267   edm::LogVerbatim("MuonSim") << " theta= " << theThetaAtEntry << " phi= " << thePhiAtEntry << " Pabs(GeV/c)  "
0268                               << thePabs << " Eloss(GeV)= " << theEnergyLoss << " Tof(ns)=  " << theTof
0269                               << " trackID= " << theTrackID << " detID= " << theDetUnitId << "\n Local:  entry "
0270                               << theEntryPoint << " exit " << theExitPoint << " delta "
0271                               << (theExitPoint - theEntryPoint) << "\n Global: entry "
0272                               << aStep->GetPreStepPoint()->GetPosition() << " exit "
0273                               << aStep->GetPostStepPoint()->GetPosition();
0274 #endif
0275 }
0276
0277 void MuonSensitiveDetector::updateHit(const G4Step* aStep) {
0278   Local3DPoint theExitPoint;
0279
0280   if (detector->isBarrel()) {
0281     theExitPoint = cmsUnits(theRotation->transformPoint(FinalStepPositionVsParent(aStep, 1), aStep));
0282   } else if (detector->isEndcap()) {
0283     // save local z at current level
0284     theExitPoint = theRotation->transformPoint(FinalStepPosition(aStep, LocalCoordinates), aStep);
0285     float zexit = theExitPoint.z();
0286     Local3DPoint tempExitPoint = theRotation->transformPoint(FinalStepPositionVsParent(aStep, 4), aStep);
0287     theExitPoint = cmsUnits(Local3DPoint(tempExitPoint.x(), tempExitPoint.y(), zexit));
0288   } else {
0289     theExitPoint = cmsUnits(theRotation->transformPoint(FinalStepPosition(aStep, LocalCoordinates), aStep));
0290   }
0291
0292   float theEnergyLoss = aStep->GetTotalEnergyDeposit() / CLHEP::GeV;
0293
0294   theHit->updateExitPoint(theExitPoint);
0295   theHit->addEnergyLoss(theEnergyLoss);
0296
0297 #ifdef EDM_ML_DEBUG
0298   edm::LogVerbatim("MuonSim") << "=== NEW Update muon hit for " << GetName() << " Edep(GeV)= " << theEnergyLoss << " "
0299                               << thePV->GetLogicalVolume()->GetName();
0300   const G4VProcess* p = aStep->GetPostStepPoint()->GetProcessDefinedStep();
0301   const G4VProcess* p2 = aStep->GetPreStepPoint()->GetProcessDefinedStep();
0302   G4String sss = "";
0303   if (p)
0304     sss += " POST PROCESS: " + p->GetProcessName();
0305   if (p2)
0306     sss += ";  PRE  PROCESS: " + p2->GetProcessName();
0307   if (!sss.empty())
0308     edm::LogVerbatim("MuonSim") << sss;
0309   edm::LogVerbatim("MuonSim") << " delEloss(GeV)= " << theEnergyLoss
0310                               << " Tof(ns)=  " << aStep->GetPreStepPoint()->GetGlobalTime() / CLHEP::nanosecond
0311                               << " trackID= " << theTrackID << " detID= " << theDetUnitId << " exit " << theExitPoint;
0312 #endif
0313 }
0314
0315 void MuonSensitiveDetector::saveHit() {
0316   if (theHit) {
0317     if (acceptHit(theHit->detUnitId())) {
0318       if (printHits_) {
0319         thePrinter->startNewSimHit(detector->name());
0320         thePrinter->printId(theHit->detUnitId());
0321         thePrinter->printLocal(theHit->entryPoint(), theHit->exitPoint());
0322       }
0323       // hit is included into hit collection
0324       slaveMuon->processHits(*theHit);
0325     }
0326     delete theHit;
0327     theHit = nullptr;
0328   }
0329 }
0330
0331 void MuonSensitiveDetector::EndOfEvent(G4HCofThisEvent*) { saveHit(); }
0332
0333 void MuonSensitiveDetector::fillHits(edm::PSimHitContainer& cc, const std::string& hname) {
0334   if (slaveMuon->name() == hname) {
0335     cc = slaveMuon->hits();
0336   }
0337 }
0338
0339 Local3DPoint MuonSensitiveDetector::InitialStepPositionVsParent(const G4Step* currentStep, G4int levelsUp) {
0340   const G4StepPoint* preStepPoint = currentStep->GetPreStepPoint();
0341   const G4ThreeVector& globalCoordinates = preStepPoint->GetPosition();
0342
0343   const G4TouchableHistory* theTouchable = (const G4TouchableHistory*)(preStepPoint->GetTouchable());
0344
0345   G4int depth = theTouchable->GetHistory()->GetDepth();
0346   G4ThreeVector localCoordinates =
0347       theTouchable->GetHistory()->GetTransform(depth - levelsUp).TransformPoint(globalCoordinates);
0348
0349   return ConvertToLocal3DPoint(localCoordinates);
0350 }
0351
0352 Local3DPoint MuonSensitiveDetector::FinalStepPositionVsParent(const G4Step* currentStep, G4int levelsUp) {
0353   const G4StepPoint* postStepPoint = currentStep->GetPostStepPoint();
0354   const G4StepPoint* preStepPoint = currentStep->GetPreStepPoint();
0355   const G4ThreeVector& globalCoordinates = postStepPoint->GetPosition();
0356
0357   const G4TouchableHistory* theTouchable = (const G4TouchableHistory*)(preStepPoint->GetTouchable());
0358
0359   G4int depth = theTouchable->GetHistory()->GetDepth();
0360   G4ThreeVector localCoordinates =
0361       theTouchable->GetHistory()->GetTransform(depth - levelsUp).TransformPoint(globalCoordinates);
0362
0363   return ConvertToLocal3DPoint(localCoordinates);
0364 }
0365
0366 bool MuonSensitiveDetector::acceptHit(uint32_t id) {
0367   if (id == 0) {
0368 #ifdef EDM_ML_DEBUG
0369     edm::LogVerbatim("MuonSim") << "DetId " << std::hex << id << std::dec << " Flag " << false;
0370 #endif
0371     return false;
0372   }
0373   if (!removeGEMHits_.empty()) {
0374     if (DetId(id).subdetId() == MuonSubdetId::GEM) {
0375       int id0 = ((GEMDetId(id).region() - GEMDetId::minRegionId)) * 10000 + GEMDetId(id).ring() * 1000 +
0376                 GEMDetId(id).station() * 100 + GEMDetId(id).chamber();
0377       if (std::find(removeGEMHits_.begin(), removeGEMHits_.end(), id0) != removeGEMHits_.end()) {
0378 #ifdef EDM_ML_DEBUG
0379         edm::LogVerbatim("MuonSim") << "Remove Hits for DetId " << GEMDetId(id);
0380 #endif
0381         return false;
0382       }
0383     }
0384   }
0385   bool flag(true);
0386   if (haveDemo_) {
0387     int subdet = DetId(id).subdetId();
0388     if (subdet == MuonSubdetId::GEM) {
0389       if (GEMDetId(id).station() == 2)
0390         flag = demoGEM_;
0391     } else if (subdet == MuonSubdetId::RPC) {
0392       if ((RPCDetId(id).region() != 0) && (RPCDetId(id).ring() == 1) && (RPCDetId(id).station() > 2))
0393         flag = demoRPC_;
0394     }
0395   }
0396 #ifdef EDM_ML_DEBUG
0397   int subdet = DetId(id).subdetId();
0398   if (subdet == MuonSubdetId::RPC)
0399     edm::LogVerbatim("MuonSim") << "DetId " << std::hex << id << std::dec << " RPC " << RPCDetId(id) << " Flag "
0400                                 << flag;
0401   else if (subdet == MuonSubdetId::GEM)
0402     edm::LogVerbatim("MuonSim") << "DetId " << std::hex << id << std::dec << " GEM " << GEMDetId(id) << " Flag "
0403                                 << flag;
0404   else if (subdet == MuonSubdetId::ME0)
0405     edm::LogVerbatim("MuonSim") << "DetId " << std::hex << id << std::dec << " " << ME0DetId(id) << " Flag " << flag;
0406   else if (subdet == MuonSubdetId::CSC)
0407     edm::LogVerbatim("MuonSim") << "DetId " << std::hex << id << std::dec << " CSC Flag " << flag;
0408   else if (subdet == MuonSubdetId::DT)
0409     edm::LogVerbatim("MuonSim") << "DetId " << std::hex << id << std::dec << " DT Flag " << flag;
0410   else
0411     edm::LogVerbatim("MuonSim") << "DetId " << std::hex << id << std::dec << " Unknown Flag " << flag;
0412 #endif
0413   return flag;
0414 }