PPS/src/PPSPixelSD.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:     PPS
0004 // Class  :     PPSPixelSD
0005 //
0006 // Implementation:
0007 //     <Notes on implementation>
0008 //
0009 // Original Author: F.Ferro
0010 //         Created:  May 4, 2015
0011 //
0012
0013 // system include files
0014
0015 // user include files
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018
0019 #include "SimG4Core/Notification/interface/TrackInformation.h"
0020 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0021 #include "SimG4Core/Physics/interface/G4ProcessTypeEnumerator.h"
0022
0023 #include "SimDataFormats/TrackingHit/interface/UpdatablePSimHit.h"
0024 #include "SimDataFormats/SimHitMaker/interface/TrackingSlaveSD.h"
0025
0026 #include "SimG4CMS/PPS/interface/PPSPixelSD.h"
0027 #include "SimG4CMS/PPS/interface/PPSPixelOrganization.h"
0028
0029 #include "G4SDManager.hh"
0030 #include "G4Step.hh"
0031 #include "G4Track.hh"
0032 #include "G4VProcess.hh"
0033
0034 #include "G4PhysicalConstants.hh"
0035 #include <CLHEP/Units/SystemOfUnits.h>
0036
0037 PPSPixelSD::PPSPixelSD(const std::string& pname,
0038                        const SensitiveDetectorCatalog& clg,
0039                        edm::ParameterSet const& p,
0040                        SimTrackManager const* manager)
0041     : SensitiveTkDetector(pname, clg), theManager_(manager) {
0042   //Add PPS Sentitive Detector Names
0043   collectionName.insert(pname);
0044
0045   //Parameters
0046   edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("PPSPixelSD");
0047   int verbn = m_p.getUntrackedParameter<int>("Verbosity");
0048   SetVerboseLevel(verbn);
0049   slave_ = std::make_unique<TrackingSlaveSD>(pname);
0050   if (pname == "CTPPSPixelHits") {
0051     numberingScheme_ = std::make_unique<PPSPixelOrganization>();
0052   } else {
0053     edm::LogWarning("PPSSim") << "PPSPixelSD: ReadoutName " << pname << " not supported";
0054   }
0055   edm::LogVerbatim("PPSSim") << "PPSPixelSD: Instantiation completed for " << pname;
0056 }
0057
0058 PPSPixelSD::~PPSPixelSD() {}
0059
0060 bool PPSPixelSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
0061   eloss_ = aStep->GetTotalEnergyDeposit();
0062   if (eloss_ > 0.0) {
0063     eloss_ /= CLHEP::GeV;
0064     stepInfo(aStep);
0065     LogDebug("PPSSim") << "PPSPixelSD: ProcessHits: edep=" << eloss_ << " "
0066                        << theTrack_->GetDefinition()->GetParticleName();
0067     if (!hitExists())
0068       createNewHit();
0069   }
0070   return true;
0071 }
0072
0073 uint32_t PPSPixelSD::setDetUnitId(const G4Step* aStep) {
0074   return (numberingScheme_ == nullptr ? 0 : numberingScheme_->unitID(aStep));
0075 }
0076
0077 void PPSPixelSD::Initialize(G4HCofThisEvent* HCE) {
0078   LogDebug("PPSSim") << "PPSPixelSD : Initialize called for " << GetName();
0079
0080   theHC_ = new PPSPixelG4HitCollection(GetName(), collectionName[0]);
0081   G4SDManager::GetSDMpointer()->AddNewCollection(GetName(), collectionName[0]);
0082
0083   if (hcID_ < 0)
0084     hcID_ = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
0085   HCE->AddHitsCollection(hcID_, theHC_);
0086
0087   tsID_ = -2;
0088   edm::LogVerbatim("PPSSim") << "PPSPixelSD: is initialized " << GetName();
0089 }
0090
0091 void PPSPixelSD::EndOfEvent(G4HCofThisEvent*) {
0092   // here we loop over transient hits and make them persistent
0093   for (unsigned int j = 0; j < (unsigned int)theHC_->entries(); ++j) {
0094     PPSPixelG4Hit* aHit = (*theHC_)[j];
0095 #ifdef debug
0096     LogDebug("PPSSim") << "HIT NUMERO " << j << "unit ID = " << aHit->unitID() << "\n"
0097                        << "               "
0098                        << "enrty z " << aHit->entry().z() << "\n"
0099                        << "               "
0100                        << "theta   " << aHit->thetaAtEntry() << "\n";
0101 #endif
0102
0103     Local3DPoint Enter(aHit->entryPoint().x(), aHit->entryPoint().y(), aHit->entryPoint().z());
0104     Local3DPoint Exit(aHit->exitPoint().x(), aHit->exitPoint().y(), aHit->exitPoint().z());
0105     slave_->processHits(PSimHit(Enter,
0106                                 Exit,
0107                                 aHit->p(),
0108                                 aHit->tof(),
0109                                 aHit->energyLoss(),
0110                                 aHit->particleType(),
0111                                 aHit->unitID(),
0112                                 aHit->trackID(),
0113                                 aHit->thetaAtEntry(),
0114                                 aHit->phiAtEntry()));
0115   }
0116 }
0117
0118 void PPSPixelSD::PrintAll() {
0119   LogDebug("PPSSim") << "PPSPixelSD: Collection " << theHC_->GetName();
0120   theHC_->PrintAllHits();
0121 }
0122
0123 void PPSPixelSD::fillHits(edm::PSimHitContainer& c, const std::string& n) {
0124   if (slave_->name() == n) {
0125     c = slave_->hits();
0126   }
0127 }
0128
0129 void PPSPixelSD::update(const BeginOfEvent* i) {
0130   LogDebug("PPSSim") << "PPSPixelSD: dispatched BeginOfEvent for " << GetName();
0131   clearHits();
0132   eventno_ = (*i)()->GetEventID();
0133 }
0134
0135 void PPSPixelSD::update(const ::EndOfEvent*) {}
0136
0137 void PPSPixelSD::clearHits() { slave_->Initialize(); }
0138
0139 G4ThreeVector PPSPixelSD::setToLocal(const G4ThreeVector& global) {
0140   G4ThreeVector localPoint;
0141   const G4VTouchable* touch = preStepPoint_->GetTouchable();
0142   localPoint = touch->GetHistory()->GetTopTransform().TransformPoint(global);
0143   return localPoint;
0144 }
0145
0146 void PPSPixelSD::stepInfo(const G4Step* aStep) {
0147   preStepPoint_ = aStep->GetPreStepPoint();
0148   postStepPoint_ = aStep->GetPostStepPoint();
0149   theTrack_ = aStep->GetTrack();
0150
0151   hitPoint_ = preStepPoint_->GetPosition();
0152   exitPoint_ = postStepPoint_->GetPosition();
0153   currentPV_ = preStepPoint_->GetPhysicalVolume();
0154   theLocalEntryPoint_ = setToLocal(hitPoint_);
0155   theLocalExitPoint_ = setToLocal(exitPoint_);
0156
0157   tSlice_ = postStepPoint_->GetGlobalTime() / CLHEP::nanosecond;
0158   tSliceID_ = (int)tSlice_;
0159   unitID_ = setDetUnitId(aStep);
0160 #ifdef debug
0161   LogDebug("PPSSim") << "UNIT ID " << unitID_;
0162 #endif
0163   primaryID_ = theTrack_->GetTrackID();
0164   parentID_ = theTrack_->GetParentID();
0165
0166   incidentEnergy_ = theTrack_->GetTotalEnergy() / CLHEP::GeV;
0167
0168   pabs_ = aStep->GetPreStepPoint()->GetMomentum().mag() / CLHEP::GeV;
0169   tof_ = aStep->GetPostStepPoint()->GetGlobalTime() / CLHEP::nanosecond;
0170
0171   particleType_ = theTrack_->GetDefinition()->GetPDGEncoding();
0172
0173   thetaAtEntry_ = aStep->GetPreStepPoint()->GetPosition().theta();
0174   phiAtEntry_ = aStep->GetPreStepPoint()->GetPosition().phi();
0175
0176   vx_ = theTrack_->GetVertexPosition().x();
0177   vy_ = theTrack_->GetVertexPosition().y();
0178   vz_ = theTrack_->GetVertexPosition().z();
0179 }
0180
0181 bool PPSPixelSD::hitExists() {
0182   // Update if in the same detector, time-slice and for same track
0183   if (tSliceID_ == tsID_ && unitID_ == previousUnitID_) {
0184     updateHit();
0185     return true;
0186   }
0187
0188   // look in the HitContainer whether a hit with the same unitID_,
0189   // tSliceID_ already exists: secondary energy deposition
0190   // is added to the primary hit
0191   int nhits = theHC_->entries();
0192   for (int j = 0; j < nhits; ++j) {
0193     PPSPixelG4Hit* aPreviousHit = (*theHC_)[j];
0194     if (aPreviousHit->timeSliceID() == tSliceID_ && aPreviousHit->unitID() == unitID_) {
0195       currentHit_ = aPreviousHit;
0196       updateHit();
0197       return true;
0198     }
0199   }
0200   return false;
0201 }
0202
0203 void PPSPixelSD::createNewHit() {
0204 #ifdef debug
0205   LogDebug("PPSSim") << "PPSPixelSD::CreateNewHit for"
0206                      << " PV " << currentPV_->GetName() << " PVid = " << currentPV_->GetCopyNo()
0207                      << " MVid = " << currentPV_->GetMother()->GetCopyNo() << " Unit " << unitID_ << "\n"
0208                      << " primary " << primaryID_ << " time slice " << tSliceID_ << " For Track  "
0209                      << " which is a " << theTrack_->GetDefinition()->GetParticleName();
0210
0211   if (parentID == 0) {
0212     LogDebug("PPSSim") << "PPSPixelSD: primary of energy " << incidentEnergy_;
0213   } else {
0214     LogDebug("PPSSim") << " daughter of track: " << parentID;
0215   }
0216
0217   if (theTrack_->GetCreatorProcess() != nullptr)
0218     LogDebug("PPSSim") << theTrack_->GetCreatorProcess()->GetProcessName();
0219   else
0220     LogDebug("PPSSim") << "NO process";
0221 #endif
0222
0223   currentHit_ = new PPSPixelG4Hit;
0224   currentHit_->setTrackID(primaryID_);
0225   currentHit_->setTimeSlice(tSlice_);
0226   currentHit_->setUnitID(unitID_);
0227   currentHit_->setIncidentEnergy(incidentEnergy_);
0228
0229   currentHit_->setP(pabs_);
0230   currentHit_->setTof(tof_);
0231   currentHit_->setParticleType(particleType_);
0232   currentHit_->setThetaAtEntry(thetaAtEntry_);
0233   currentHit_->setPhiAtEntry(phiAtEntry_);
0234
0235   G4ThreeVector pos = 0.5 * (hitPoint_ + exitPoint_);
0236   currentHit_->setMeanPosition(pos);
0237   currentHit_->setEntryPoint(theLocalEntryPoint_);
0238   currentHit_->setExitPoint(theLocalExitPoint_);
0239
0240   currentHit_->setParentId(parentID_);
0241   currentHit_->setVx(vx_);
0242   currentHit_->setVy(vy_);
0243   currentHit_->setVz(vz_);
0244
0245   updateHit();
0246   storeHit(currentHit_);
0247 }
0248
0249 void PPSPixelSD::updateHit() {
0250 #ifdef debug
0251   LogDebug("PPSSim") << "G4PPSPixelSD updateHit: add eloss=" << eloss_ << "\nCurrentHit=" << currentHit_
0252                      << ", PostStepPoint=" << postStepPoint_->GetPosition();
0253 #endif
0254   currentHit_->setEnergyLoss(eloss_);
0255   // buffer for next steps:
0256   tsID_ = tSliceID_;
0257   previousUnitID_ = unitID_;
0258 }
0259
0260 void PPSPixelSD::storeHit(PPSPixelG4Hit* hit) {
0261   if (hit == nullptr) {
0262     edm::LogWarning("PPSSim") << "PPSPixelSD: hit to be stored is NULL !!";
0263     return;
0264   }
0265
0266   theHC_->insert(hit);
0267 }