MuonSensitiveDetector.cc

Go to the documentation of this file.

#include "SimG4CMS/Muon/interface/MuonSensitiveDetector.h"
#include "SimG4CMS/Muon/interface/MuonSlaveSD.h"
#include "SimG4CMS/Muon//interface/MuonEndcapFrameRotation.h"
#include "SimG4CMS/Muon/interface/MuonRPCFrameRotation.h"
#include "SimG4CMS/Muon/interface/MuonGEMFrameRotation.h"
#include "SimG4CMS/Muon/interface/MuonME0FrameRotation.h"
#include "Geometry/MuonNumbering/interface/MuonSubDetector.h"
#include "Geometry/MuonNumbering/interface/MuonDDDConstants.h"

#include "SimG4CMS/Muon/interface/SimHitPrinter.h"
#include "SimDataFormats/TrackingHit/interface/UpdatablePSimHit.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"

#include "SimG4CMS/Muon/interface/MuonG4Numbering.h"
#include "Geometry/MuonNumbering/interface/MuonBaseNumber.h"
#include "Geometry/MuonNumbering/interface/MuonSimHitNumberingScheme.h"

#include "SimG4Core/Notification/interface/TrackInformation.h"
#include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
#include "SimG4Core/Physics/interface/G4ProcessTypeEnumerator.h"

#include "G4VProcess.hh"
#include "G4EventManager.hh"
#include "G4SystemOfUnits.hh"
#include "G4Step.hh"
#include "G4StepPoint.hh"
#include "G4Track.hh"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include <iostream>

//#define DebugLog

MuonSensitiveDetector::MuonSensitiveDetector(const std::string& name, 
                         const DDCompactView & cpv,
                         const SensitiveDetectorCatalog & clg,
                         edm::ParameterSet const & p,
                         const SimTrackManager* manager) 
  : SensitiveTkDetector(name, cpv, clg, p),
    thePV(nullptr), theHit(nullptr), theDetUnitId(0), newDetUnitId(0), 
    theTrackID(0), theManager(manager)
{
  edm::ParameterSet m_MuonSD = p.getParameter<edm::ParameterSet>("MuonSD");
  ePersistentCutGeV = m_MuonSD.getParameter<double>("EnergyThresholdForPersistency")/CLHEP::GeV;//Default 1. GeV
  allMuonsPersistent = m_MuonSD.getParameter<bool>("AllMuonsPersistent");
  printHits  = m_MuonSD.getParameter<bool>("PrintHits");
  
  //
  // Here simply create 1 MuonSlaveSD for the moment
  //  
  LogDebug("MuonSimDebug") << "create MuonSubDetector "<<name;
  detector = new MuonSubDetector(name);

  //The constants take time to calculate and are needed by many helpers
  MuonDDDConstants constants(cpv);
  G4String sdet = "unknown";
  if (detector->isEndcap()) {
    theRotation=new MuonEndcapFrameRotation();
    sdet = "Endcap";
  } else if (detector->isRPC()) {
    theRotation=new MuonRPCFrameRotation( constants );
    sdet = "RPC";
  } else if (detector->isGEM()) {
    theRotation=new MuonGEMFrameRotation( constants );
    sdet = "GEM";
  } else if (detector->isME0()) {
    theRotation=new MuonME0FrameRotation( constants );
    sdet = "ME0";
  } else {
    theRotation = new MuonFrameRotation();
  }
  slaveMuon  = new MuonSlaveSD(detector,theManager);
  numbering  = new MuonSimHitNumberingScheme(detector, constants);
  g4numbering = new MuonG4Numbering(constants);
  
  if (printHits) {
    thePrinter = new SimHitPrinter("HitPositionOSCAR.dat");
  }

  edm::LogVerbatim("MuonSensitiveDetector") 
    << " of type " << sdet << " <" << GetName() 
    << "> EnergyThresholdForPersistency(GeV) " << ePersistentCutGeV/CLHEP::GeV 
    << " allMuonsPersistent: " << allMuonsPersistent;
    
  theG4ProcessTypeEnumerator = new G4ProcessTypeEnumerator;
}

MuonSensitiveDetector::~MuonSensitiveDetector() { 
  delete g4numbering;
  delete numbering;
  delete slaveMuon;
  delete theRotation;
  delete detector;
  delete theG4ProcessTypeEnumerator;
}

void MuonSensitiveDetector::update(const BeginOfEvent * i){
  clearHits();
  //----- Initialize variables to check if two steps belong to same hit
  thePV = nullptr;
  theDetUnitId = 0;
  theTrackID = 0;
}

void MuonSensitiveDetector::clearHits()
{
  LogDebug("MuonSimDebug") << "MuonSensitiveDetector::clearHits";
  slaveMuon->Initialize();
}

bool MuonSensitiveDetector::ProcessHits(G4Step * aStep, G4TouchableHistory * ROhist)
{
  LogDebug("MuonSimDebug") <<" MuonSensitiveDetector::ProcessHits "
               <<InitialStepPosition(aStep,WorldCoordinates);

  if (aStep->GetTotalEnergyDeposit()>0.){
    newDetUnitId = setDetUnitId(aStep);

    // do not count neutrals that are killed by User Limits MinEKine
    //---VI: this is incorrect, neutral particle, like neutron may have local
    //       energy deposit, which potentially may make a hit
    if( aStep->GetTrack()->GetDynamicParticle()->GetCharge() != 0 ){
  
      if (newHit(aStep)) {
    saveHit();
    createHit(aStep);
      } else {
    updateHit(aStep);
      } 
    } else {
      thePV = aStep->GetPreStepPoint()->GetPhysicalVolume();
      theTrackID = aStep->GetTrack()->GetTrackID();
      theDetUnitId = newDetUnitId;
    }
  }
  return true;
}

uint32_t MuonSensitiveDetector::setDetUnitId(const G4Step * aStep)
{ 
  MuonBaseNumber num = g4numbering->PhysicalVolumeToBaseNumber(aStep);

#ifdef DebugLog
  std::stringstream MuonBaseNumber; 
  MuonBaseNumber << "MuonNumbering :: number of levels = "<<num.getLevels()<<std::endl;
  MuonBaseNumber << "Level \t SuperNo \t BaseNo"<<std::endl;
  for (int level=1;level<=num.getLevels();level++) {
    MuonBaseNumber << level << " \t " << num.getSuperNo(level)
          << " \t " << num.getBaseNo(level) << std::endl;
  }
  std::string MuonBaseNumbr = MuonBaseNumber.str();

  LogDebug("MuonSimDebug") <<"MuonSensitiveDetector::setDetUnitId :: "<<MuonBaseNumbr;
  LogDebug("MuonSimDebug") <<"MuonSensitiveDetector::setDetUnitId :: MuonDetUnitId = "
               <<(numbering->baseNumberToUnitNumber(num));
#endif
  return numbering->baseNumberToUnitNumber(num);
}

bool MuonSensitiveDetector::newHit(const G4Step * aStep){

  return (!theHit || (aStep->GetTrack()->GetTrackID() != theTrackID) 
      || (aStep->GetPreStepPoint()->GetPhysicalVolume() != thePV)
      || newDetUnitId != theDetUnitId);
}

void MuonSensitiveDetector::createHit(const G4Step * aStep){

  Local3DPoint theEntryPoint;
  Local3DPoint theExitPoint;

  if (detector->isBarrel()) {
    // 1 levels up
    theEntryPoint = cmsUnits(theRotation->transformPoint(InitialStepPositionVsParent(aStep,1),aStep));
    theExitPoint  = cmsUnits(theRotation->transformPoint(FinalStepPositionVsParent(aStep,1),aStep));
  } else if (detector->isEndcap()) {
    // save local z at current level
    theEntryPoint = theRotation->transformPoint(InitialStepPosition(aStep,LocalCoordinates),aStep);
    theExitPoint  = theRotation->transformPoint(FinalStepPosition(aStep,LocalCoordinates),aStep);
    float zentry  = theEntryPoint.z();
    float zexit   = theExitPoint.z();
    // 4 levels up
    Local3DPoint tempEntry = theRotation->transformPoint(InitialStepPositionVsParent(aStep,4),aStep);
    Local3DPoint tempExit  = theRotation->transformPoint(FinalStepPositionVsParent(aStep,4),aStep);
    // reset local z from z wrt deep-parent volume to z wrt low-level volume
    theEntryPoint = cmsUnits(Local3DPoint( tempEntry.x(), tempEntry.y(), zentry ));
    theExitPoint  = cmsUnits(Local3DPoint( tempExit.x(),  tempExit.y(),  zexit ));
  } else {
    theEntryPoint = cmsUnits(theRotation->transformPoint(InitialStepPosition(aStep,LocalCoordinates),aStep));
    theExitPoint  = cmsUnits(theRotation->transformPoint(FinalStepPosition(aStep,LocalCoordinates),aStep)); 
  }

  const G4Track* theTrack = aStep->GetTrack(); 
  const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();

  float thePabs       = preStepPoint->GetMomentum().mag()/CLHEP::GeV;
  float theTof        = preStepPoint->GetGlobalTime()/CLHEP::nanosecond;
  float theEnergyLoss = aStep->GetTotalEnergyDeposit()/CLHEP::GeV;
  int theParticleType = G4TrackToParticleID::particleID(theTrack);

  theDetUnitId = newDetUnitId;
  thePV = preStepPoint->GetPhysicalVolume();
  theTrackID = theTrack->GetTrackID();

  // convert momentum direction it to local frame
  const G4ThreeVector& gmd  = preStepPoint->GetMomentumDirection();
  G4ThreeVector lmd = ((G4TouchableHistory *)(preStepPoint->GetTouchable()))->GetHistory()
      ->GetTopTransform().TransformAxis(gmd);
  Local3DPoint lnmd = ConvertToLocal3DPoint(lmd);
  lnmd = theRotation->transformPoint(lnmd, aStep);
  float theThetaAtEntry = lnmd.theta();
  float thePhiAtEntry = lnmd.phi();

  theHit = new UpdatablePSimHit(theEntryPoint,theExitPoint,thePabs,theTof,
                theEnergyLoss,theParticleType,theDetUnitId,
                theTrackID,theThetaAtEntry,thePhiAtEntry,
                theG4ProcessTypeEnumerator->processId(theTrack->GetCreatorProcess()));

  // Make track persistent
  int thePID = std::abs(theTrack->GetDefinition()->GetPDGEncoding());
  //---VI - in parameters cut in energy is declared but applied to momentum
  if(thePabs > ePersistentCutGeV || ( thePID == 13 && allMuonsPersistent ) ){
    TrackInformation* info = cmsTrackInformation(theTrack);
    info->storeTrack(true);
  }

#ifdef DebugLog
  edm::LogVerbatim("MuonSimDebug") <<"=== NEW Muon hit for "<< GetName() 
                   << " Edep(GeV)= " << theEnergyLoss
                   <<" " <<thePV->GetLogicalVolume()->GetName();
  const G4VProcess* p = aStep->GetPostStepPoint()->GetProcessDefinedStep();
  const G4VProcess* p2 = preStepPoint->GetProcessDefinedStep();
  G4String sss = "";
  if (p) sss += " POST PROCESS: " + p->GetProcessName();
  if (p2)sss += ";  PRE  PROCESS: " + p2->GetProcessName();
  if("" != sss) edm::LogVerbatim("MuonSimDebug") << sss;
  edm::LogVerbatim("MuonSimDebug") << " theta= " << theThetaAtEntry
                   << " phi= " << thePhiAtEntry
                   << " Pabs(GeV/c)  " << thePabs
                   << " Eloss(GeV)= " << theEnergyLoss
                   << " Tof(ns)=  " << theTof
                   << " trackID= " << theTrackID
                   << " detID= " << theDetUnitId
                   << "\n Local:  entry " << theEntryPoint 
                   << " exit " << theExitPoint
                   << " delta " << (theExitPoint-theEntryPoint)
                   << "\n Global: entry " << preStepPoint->GetPosition() 
                   << " exit " << aStep->GetPostStepPoint()->GetPosition();
#endif
}

void MuonSensitiveDetector::updateHit(const G4Step * aStep){

  Local3DPoint theExitPoint;

  if (detector->isBarrel()) {
    theExitPoint = cmsUnits(theRotation->transformPoint(FinalStepPositionVsParent(aStep,1),aStep));  
  } else if (detector->isEndcap()) {
    // save local z at current level
    theExitPoint = theRotation->transformPoint(FinalStepPosition(aStep,LocalCoordinates),aStep);
    float zexit  = theExitPoint.z();
    Local3DPoint tempExitPoint = theRotation->transformPoint(FinalStepPositionVsParent(aStep,4),aStep);
    theExitPoint = cmsUnits(Local3DPoint( tempExitPoint.x(), tempExitPoint.y(), zexit));
  } else {
    theExitPoint = cmsUnits(theRotation->transformPoint(FinalStepPosition(aStep,LocalCoordinates),aStep)); 
  }

  float theEnergyLoss = aStep->GetTotalEnergyDeposit()/CLHEP::GeV;  

  theHit->updateExitPoint(theExitPoint);
  theHit->addEnergyLoss(theEnergyLoss);

#ifdef DebugLog
  edm::LogVerbatim("MuonSimDebug") <<"=== NEW Update muon hit for "<< GetName() 
                   << " Edep(GeV)= " << theEnergyLoss
                   <<" " <<thePV->GetLogicalVolume()->GetName();
  const G4VProcess* p = aStep->GetPostStepPoint()->GetProcessDefinedStep();
  const G4VProcess* p2 = preStepPoint->GetProcessDefinedStep();
  G4String sss = "";
  if (p) sss += " POST PROCESS: " + p->GetProcessName();
  if (p2)sss += ";  PRE  PROCESS: " + p2->GetProcessName();
  if("" != sss) edm::LogVerbatim("MuonSimDebug") << sss;
  edm::LogVerbatim("MuonSimDebug") << " delEloss(GeV)= " << theEnergyLoss
                   << " Tof(ns)=  " << theTof
                   << " trackID= " << theTrackID
                   << " detID= " << theDetUnitId
                   << " exit " << theExitPoint;
#endif
}

void MuonSensitiveDetector::saveHit(){
  if (theHit) {
    if (printHits) {
      thePrinter->startNewSimHit(detector->name());
      thePrinter->printId(theHit->detUnitId());
      thePrinter->printLocal(theHit->entryPoint(),theHit->exitPoint());
    }
    // hit is included into hit collection
    slaveMuon->processHits(*theHit);
    delete theHit;
    theHit = nullptr;
  }
}

void MuonSensitiveDetector::EndOfEvent(G4HCofThisEvent*)
{
  saveHit();
}

void MuonSensitiveDetector::fillHits(edm::PSimHitContainer& cc, const std::string& hname){
  if (slaveMuon->name() == hname) { cc=slaveMuon->hits(); }
}

Local3DPoint MuonSensitiveDetector::InitialStepPositionVsParent(const G4Step * currentStep, G4int levelsUp) {
  
  const G4StepPoint * preStepPoint = currentStep->GetPreStepPoint();
  const G4ThreeVector& globalCoordinates = preStepPoint->GetPosition();
  
  const G4TouchableHistory * theTouchable=(const G4TouchableHistory *)
    (preStepPoint->GetTouchable());

  G4int depth = theTouchable->GetHistory()->GetDepth();
  G4ThreeVector localCoordinates = theTouchable->GetHistory()
    ->GetTransform(depth-levelsUp).TransformPoint(globalCoordinates);
  
  return ConvertToLocal3DPoint(localCoordinates); 
}
 
Local3DPoint MuonSensitiveDetector::FinalStepPositionVsParent(const G4Step* currentStep, G4int levelsUp) {
  
  const G4StepPoint * postStepPoint = currentStep->GetPostStepPoint();
  const G4StepPoint * preStepPoint  = currentStep->GetPreStepPoint();
  const G4ThreeVector& globalCoordinates = postStepPoint->GetPosition();
    
  const G4TouchableHistory * theTouchable = (const G4TouchableHistory *)
    (preStepPoint->GetTouchable());

  G4int depth = theTouchable->GetHistory()->GetDepth();
  G4ThreeVector localCoordinates = theTouchable->GetHistory()
    ->GetTransform(depth-levelsUp).TransformPoint(globalCoordinates);

  return ConvertToLocal3DPoint(localCoordinates); 
}