PPS/src/PPSDiamondSD.cc

0001 ///////////////////////////////////////////////////////////////////////////////
0002 //Author: Seyed Mohsen Etesami
0003 // setesami@cern.ch
0004 // 2016 Nov
0005 ///////////////////////////////////////////////////////////////////////////////
0006 #include "SimG4CMS/PPS/interface/PPSDiamondSD.h"
0007 #include "SimG4CMS/PPS/interface/PPSDiamondOrganization.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "SimG4Core/Notification/interface/TrackInformation.h"
0011 #include "SimG4Core/Geometry/interface/SensitiveDetectorCatalog.h"
0012 #include "SimG4Core/Physics/interface/G4ProcessTypeEnumerator.h"
0013 #include "SimG4Core/SensitiveDetector/interface/SensitiveTkDetector.h"
0014 #include "SimDataFormats/TrackingHit/interface/UpdatablePSimHit.h"
0015 #include "SimDataFormats/SimHitMaker/interface/TrackingSlaveSD.h"
0016 #include "G4SDManager.hh"
0017 #include "G4Step.hh"
0018 #include "G4Track.hh"
0019 #include "G4VProcess.hh"
0020 #include "G4ParticleDefinition.hh"
0021 #include "G4ParticleTypes.hh"
0022 #include "SimG4Core/Notification/interface/TrackInformation.h"
0023 #include "G4PhysicalConstants.hh"
0024 #include <CLHEP/Units/SystemOfUnits.h>
0025 #include <iostream>
0026 #include <vector>
0027 #include <string>
0028
0029 PPSDiamondSD::PPSDiamondSD(const std::string& pname,
0030                            const SensitiveDetectorCatalog& clg,
0031                            edm::ParameterSet const& p,
0032                            const SimTrackManager* manager)
0033     : SensitiveTkDetector(pname, clg) {
0034   collectionName.insert(pname);
0035   edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("PPSDiamondSD");
0036   verbosity_ = m_Anal.getParameter<int>("Verbosity");
0037
0038   LogDebug("PPSSimDiamond") << "*******************************************************\n"
0039                             << "*                                                     *\n"
0040                             << "* Constructing a PPSDiamondSD  with name " << pname << "\n"
0041                             << "*                                                     *\n"
0042                             << "*******************************************************"
0043                             << "\n";
0044
0045   slave_ = std::make_unique<TrackingSlaveSD>(pname);
0046
0047   if (pname == "CTPPSTimingHits") {
0048     numberingScheme_ = std::make_unique<PPSDiamondOrganization>();
0049     edm::LogVerbatim("PPSSimDiamond") << "Find CTPPSTimingHits as name";
0050   } else {
0051     edm::LogError("PPSSimDiamond") << "PPSDiamondSD: ReadoutName " << pname << " not supported";
0052   }
0053
0054   edm::LogVerbatim("PPSSimDiamond") << "PPSDiamondSD: Instantiation completed for " << pname;
0055 }
0056
0057 PPSDiamondSD::~PPSDiamondSD() {}
0058
0059 void PPSDiamondSD::Initialize(G4HCofThisEvent* HCE) {
0060   LogDebug("PPSSimDiamond") << "PPSDiamondSD : Initialize called for " << GetName();
0061
0062   theHC_ = new PPSDiamondG4HitCollection(GetName(), collectionName[0]);
0063   G4SDManager::GetSDMpointer()->AddNewCollection(GetName(), collectionName[0]);
0064   if (hcID_ < 0)
0065     hcID_ = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
0066   HCE->AddHitsCollection(hcID_, theHC_);
0067 }
0068
0069 void PPSDiamondSD::printHitInfo() {
0070   LogDebug("PPSSimDiamond") << theTrack_->GetDefinition()->GetParticleName() << " PPS_Timing_SD CreateNewHit for"
0071                             << " PV " << currentPV_->GetName() << " PVid = " << currentPV_->GetCopyNo() << " Unit "
0072                             << unitID_ << "\n";
0073   LogDebug("PPSSimDiamond") << " primary " << primaryID_ << " time slice " << tSliceID_ << " of energy "
0074                             << theTrack_->GetTotalEnergy() << " Eloss=" << eloss_ << " position pre: " << hitPoint_
0075                             << " post: " << exitPoint_;
0076   LogDebug("PPSSimDiamond") << " positions "
0077                             << "(" << postStepPoint_->GetPosition().x() << "," << postStepPoint_->GetPosition().y()
0078                             << "," << postStepPoint_->GetPosition().z() << ")"
0079                             << " For Track  " << theTrack_->GetTrackID() << " which is a "
0080                             << theTrack_->GetDefinition()->GetParticleName() << " ParentID is "
0081                             << theTrack_->GetParentID() << "\n";
0082
0083   if (parentID_ == 0) {
0084     LogDebug("PPSSimDiamond") << " primary particle ";
0085   } else {
0086     LogDebug("PPSSimDiamond") << " daughter of part. " << parentID_;
0087   }
0088
0089   LogDebug("PPSSimDiamond") << " and created by ";
0090
0091   if (theTrack_->GetCreatorProcess() != nullptr)
0092     LogDebug("PPSSimDiamond") << theTrack_->GetCreatorProcess()->GetProcessName();
0093   else
0094     LogDebug("PPSSimDiamond") << "NO process";
0095 }
0096
0097 bool PPSDiamondSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
0098   eloss_ = aStep->GetTotalEnergyDeposit();
0099   if (eloss_ > 0.0) {
0100     eloss_ /= CLHEP::GeV;
0101     stepInfo(aStep);
0102     if (theTrack_->GetDefinition()->GetPDGEncoding() == 2212) {
0103       LogDebug("PPSSimDiamond") << "PPSSimDiamond : eloss=" << eloss_;
0104       importInfoToHit();  //in addition to import info to hit it STORE hit as well
0105       LogDebug("PPSSimDiamond") << " information imported to the hit ";
0106     }
0107   }
0108   return true;
0109 }
0110
0111 void PPSDiamondSD::stepInfo(const G4Step* aStep) {
0112   theTrack_ = aStep->GetTrack();
0113   preStepPoint_ = aStep->GetPreStepPoint();
0114   postStepPoint_ = aStep->GetPostStepPoint();
0115   hitPoint_ = preStepPoint_->GetPosition();
0116   exitPoint_ = postStepPoint_->GetPosition();
0117   currentPV_ = preStepPoint_->GetPhysicalVolume();
0118   theLocalEntryPoint_ = setToLocal(hitPoint_);
0119   theLocalExitPoint_ = setToLocal(exitPoint_);
0120   tof_ = preStepPoint_->GetGlobalTime() / CLHEP::nanosecond;
0121   incidentEnergy_ = preStepPoint_->GetTotalEnergy() / CLHEP::eV;
0122   tSlice_ = postStepPoint_->GetGlobalTime() / CLHEP::nanosecond;
0123   tSliceID_ = (int)tSlice_;
0124   unitID_ = setDetUnitId(aStep);
0125
0126   if (verbosity_)
0127     LogDebug("PPSSimDiamond") << "UNIT " << unitID_;
0128
0129   primaryID_ = theTrack_->GetTrackID();
0130   parentID_ = theTrack_->GetParentID();
0131
0132   incidentEnergy_ = theTrack_->GetTotalEnergy() / CLHEP::GeV;
0133
0134   pabs_ = preStepPoint_->GetMomentum().mag() / CLHEP::GeV;
0135   thePx_ = preStepPoint_->GetMomentum().x() / CLHEP::GeV;
0136   thePy_ = preStepPoint_->GetMomentum().y() / CLHEP::GeV;
0137   thePz_ = preStepPoint_->GetMomentum().z() / CLHEP::GeV;
0138   particleType_ = theTrack_->GetDefinition()->GetPDGEncoding();
0139
0140   //corrected phi and theta treatment
0141   G4ThreeVector gmd = preStepPoint_->GetMomentumDirection();
0142   // convert it to local frame
0143   G4ThreeVector lmd =
0144       ((G4TouchableHistory*)(preStepPoint_->GetTouchable()))->GetHistory()->GetTopTransform().TransformAxis(gmd);
0145   Local3DPoint lnmd = ConvertToLocal3DPoint(lmd);
0146   thetaAtEntry_ = lnmd.theta();
0147   phiAtEntry_ = lnmd.phi();
0148
0149   vx_ = theTrack_->GetVertexPosition().x() / CLHEP::mm;
0150   vy_ = theTrack_->GetVertexPosition().y() / CLHEP::mm;
0151   vz_ = theTrack_->GetVertexPosition().z() / CLHEP::mm;
0152 }
0153
0154 uint32_t PPSDiamondSD::setDetUnitId(const G4Step* aStep) {
0155   return (numberingScheme_ == nullptr ? 0 : numberingScheme_->unitID(aStep));
0156 }
0157
0158 void PPSDiamondSD::storeHit(PPSDiamondG4Hit* hit) {
0159   if (hit == nullptr) {
0160     if (verbosity_)
0161       LogDebug("PPSSimDiamond") << "PPSDiamond: hit to be stored is NULL !!";
0162     return;
0163   }
0164
0165   theHC_->insert(hit);
0166 }
0167
0168 void PPSDiamondSD::importInfoToHit() {
0169   currentHit_ = new PPSDiamondG4Hit;
0170   currentHit_->setTrackID(primaryID_);
0171   currentHit_->setTimeSlice(tSlice_);
0172   currentHit_->setUnitID(unitID_);
0173   currentHit_->setIncidentEnergy(incidentEnergy_);
0174   currentHit_->setP(pabs_);
0175   currentHit_->setTof(tof_);
0176   currentHit_->setEnergyLoss(eloss_);
0177   currentHit_->setParticleType(particleType_);
0178   currentHit_->setThetaAtEntry(thetaAtEntry_);
0179   currentHit_->setPhiAtEntry(phiAtEntry_);
0180   currentHit_->setEntry(hitPoint_);
0181   currentHit_->setExit(exitPoint_);
0182   currentHit_->setLocalEntry(theLocalEntryPoint_);
0183   currentHit_->setLocalExit(theLocalExitPoint_);
0184   currentHit_->setParentId(parentID_);
0185   currentHit_->setVx(vx_);
0186   currentHit_->setVy(vy_);
0187   currentHit_->setVz(vz_);
0188   currentHit_->setPx(thePx_);
0189   currentHit_->setPy(thePy_);
0190   currentHit_->setPz(thePz_);
0191   currentHit_->setGlobalTimehit(tof_);
0192
0193   storeHit(currentHit_);
0194   LogDebug("PPSSimDiamond") << "STORED HIT IN: " << unitID_ << "\n";
0195 }
0196
0197 G4ThreeVector PPSDiamondSD::setToLocal(const G4ThreeVector& global) {
0198   G4ThreeVector localPoint;
0199   const G4VTouchable* touch = preStepPoint_->GetTouchable();
0200   localPoint = touch->GetHistory()->GetTopTransform().TransformPoint(global);
0201
0202   return localPoint;
0203 }
0204
0205 void PPSDiamondSD::EndOfEvent(G4HCofThisEvent*) {
0206   // here we loop over transient hits and make them persistent
0207   for (unsigned int j = 0; j < (unsigned int)theHC_->entries(); ++j) {
0208     PPSDiamondG4Hit* aHit = (*theHC_)[j];
0209
0210     Local3DPoint entry(aHit->localEntry().x(), aHit->localEntry().y(), aHit->localEntry().z());
0211     Local3DPoint exit(aHit->localExit().x(), aHit->localExit().y(), aHit->localExit().z());
0212     slave_->processHits(PSimHit(entry,
0213                                 exit,
0214                                 aHit->p(),
0215                                 aHit->tof(),
0216                                 aHit->energyLoss(),
0217                                 aHit->particleType(),
0218                                 aHit->unitID(),
0219                                 aHit->trackID(),
0220                                 aHit->thetaAtEntry(),
0221                                 aHit->phiAtEntry()));
0222   }
0223 }
0224
0225 void PPSDiamondSD::PrintAll() {
0226   LogDebug("PPSSimDiamond") << "PPSDiamond: Collection " << theHC_->GetName() << "\n";
0227   theHC_->PrintAllHits();
0228 }
0229
0230 void PPSDiamondSD::fillHits(edm::PSimHitContainer& c, const std::string& n) {
0231   if (slave_->name() == n)
0232     c = slave_->hits();
0233 }
0234
0235 void PPSDiamondSD::update(const BeginOfEvent* i) {
0236   LogDebug("PPSSimDiamond") << " Dispatched BeginOfEvent !"
0237                             << "\n";
0238   clearHits();
0239   eventno_ = (*i)()->GetEventID();
0240 }
0241
0242 void PPSDiamondSD::update(const ::EndOfEvent*) {}
0243
0244 void PPSDiamondSD::clearHits() { slave_->Initialize(); }