Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-01-13 01:44:06 

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