File indexing completed on 2024-05-10 02:21:15
0001 #include "SimG4Core/Geometry/interface/DD4hep2DDDName.h"
0002 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0003 #include "SimG4Core/Notification/interface/TrackInformation.h"
0004 #include "SimG4Core/Notification/interface/SimTrackManager.h"
0005
0006 #include "SimG4CMS/Calo/interface/CaloTrkProcessing.h"
0007
0008 #include "FWCore/Utilities/interface/Exception.h"
0009
0010 #include "G4EventManager.hh"
0011 #include "G4LogicalVolumeStore.hh"
0012 #include "G4LogicalVolume.hh"
0013 #include "G4Step.hh"
0014 #include "G4Track.hh"
0015
0016 #include <CLHEP/Units/SystemOfUnits.h>
0017
0018 #include <sstream>
0019
0020
0021 CaloTrkProcessing::CaloTrkProcessing(const std::string& name,
0022 const CaloSimulationParameters& csps,
0023 const SensitiveDetectorCatalog& clg,
0024 bool testBeam,
0025 double eMin,
0026 bool putHistory,
0027 bool doFineCalo,
0028 double eMinFine,
0029 int addlevel,
0030 const std::vector<std::string>& fineNames,
0031 const std::vector<int>& fineLevels,
0032 const std::vector<int>& useFines,
0033 const SimTrackManager*)
0034 : SensitiveCaloDetector(name, clg),
0035 testBeam_(testBeam),
0036 eMin_(eMin),
0037 putHistory_(putHistory),
0038 doFineCalo_(doFineCalo),
0039 eMinFine_(eMinFine),
0040 addlevel_(addlevel),
0041 lastTrackID_(-1) {
0042
0043
0044 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing: Initialised with TestBeam = " << testBeam_ << " Emin = " << eMin_
0045 << " Flags " << putHistory_ << " (History), " << doFineCalo_ << " (Special Calorimeter)";
0046 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing: Have a possibility of " << fineNames.size()
0047 << " fine calorimeters of which " << useFines.size() << " are selected";
0048 for (unsigned int k = 0; k < fineNames.size(); ++k)
0049 edm::LogVerbatim("CaloSim") << "[" << k << "] " << fineNames[k] << " at " << fineLevels[k];
0050 std::ostringstream st1;
0051 for (unsigned int k = 0; k < useFines.size(); ++k)
0052 st1 << " [" << k << "] " << useFines[k] << ":" << fineNames[useFines[k]];
0053 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing used calorimeters" << st1.str();
0054
0055
0056 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing: " << csps.caloNames_.size() << " entries for caloNames:";
0057 for (unsigned int i = 0; i < csps.caloNames_.size(); i++)
0058 edm::LogVerbatim("CaloSim") << " (" << i << ") " << csps.caloNames_[i];
0059 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing: " << csps.levels_.size() << " entries for levels:";
0060 for (unsigned int i = 0; i < csps.levels_.size(); i++)
0061 edm::LogVerbatim("CaloSim") << " (" << i << ") " << (csps.levels_[i] + addlevel_);
0062 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing: " << csps.neighbours_.size() << " entries for neighbours:";
0063 for (unsigned int i = 0; i < csps.neighbours_.size(); i++)
0064 edm::LogVerbatim("CaloSim") << " (" << i << ") " << csps.neighbours_[i];
0065 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing: " << csps.insideNames_.size() << " entries for insideNames:";
0066 for (unsigned int i = 0; i < csps.insideNames_.size(); i++)
0067 edm::LogVerbatim("CaloSim") << " (" << i << ") " << csps.insideNames_[i];
0068 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing: " << csps.insideLevel_.size() << " entries for insideLevel:";
0069 for (unsigned int i = 0; i < csps.insideLevel_.size(); i++)
0070 edm::LogVerbatim("CaloSim") << " (" << i << ") " << (csps.insideLevel_[i] + addlevel_);
0071
0072 if (csps.caloNames_.size() < csps.neighbours_.size()) {
0073 edm::LogError("CaloSim") << "CaloTrkProcessing: # of Calorimeter bins " << csps.caloNames_.size()
0074 << " does not match with " << csps.neighbours_.size() << " ==> illegal ";
0075 throw cms::Exception("Unknown", "CaloTrkProcessing")
0076 << "Calorimeter array size does not match with size of neighbours\n";
0077 }
0078
0079 const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0080 std::vector<G4LogicalVolume*>::const_iterator lvcite;
0081 int istart = 0;
0082 for (unsigned int i = 0; i < csps.caloNames_.size(); i++) {
0083 G4LogicalVolume* lv = nullptr;
0084 G4String name(csps.caloNames_[i]);
0085 for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
0086 G4String namx(DD4hep2DDDName::noNameSpace(static_cast<std::string>((*lvcite)->GetName())));
0087 if (namx == name) {
0088 lv = (*lvcite);
0089 break;
0090 }
0091 }
0092 if (lv != nullptr) {
0093 CaloTrkProcessing::Detector detector;
0094 detector.name = name;
0095 detector.lv = lv;
0096 detector.level = (csps.levels_[i] + addlevel_);
0097 if (istart + csps.neighbours_[i] > static_cast<int>(csps.insideNames_.size())) {
0098 edm::LogError("CaloSim") << "CaloTrkProcessing: # of InsideNames bins " << csps.insideNames_.size()
0099 << " too few compaerd to " << istart + csps.neighbours_[i]
0100 << " requested ==> illegal ";
0101 throw cms::Exception("Unknown", "CaloTrkProcessing")
0102 << "InsideNames array size does not match with list of neighbours\n";
0103 }
0104 std::vector<std::string> inside;
0105 std::vector<G4LogicalVolume*> insideLV;
0106 std::vector<int> insideLevels;
0107 for (int k = 0; k < csps.neighbours_[i]; k++) {
0108 lv = nullptr;
0109 name = static_cast<G4String>(csps.insideNames_[istart + k]);
0110 for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
0111 G4String namx(DD4hep2DDDName::noNameSpace(static_cast<std::string>((*lvcite)->GetName())));
0112 if (namx == name) {
0113 lv = (*lvcite);
0114 break;
0115 }
0116 }
0117 inside.push_back(name);
0118 insideLV.push_back(lv);
0119 insideLevels.push_back(csps.insideLevel_[istart + k] + addlevel_);
0120 }
0121 detector.fromDets = inside;
0122 detector.fromDetL = insideLV;
0123 detector.fromLevels = insideLevels;
0124 detectors_.emplace_back(detector);
0125 }
0126 istart += csps.neighbours_[i];
0127 }
0128
0129 for (unsigned int i = 0; i < useFines.size(); i++) {
0130 G4LogicalVolume* lv = nullptr;
0131 G4String name = static_cast<G4String>(fineNames[useFines[i]]);
0132 for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
0133 G4String namx(DD4hep2DDDName::noNameSpace(static_cast<std::string>((*lvcite)->GetName())));
0134 if (namx == name) {
0135 lv = (*lvcite);
0136 break;
0137 }
0138 }
0139 if (lv != nullptr) {
0140 CaloTrkProcessing::Detector detector;
0141 detector.name = name;
0142 detector.lv = lv;
0143 detector.level = fineLevels[useFines[i]];
0144 detector.fromDets.clear();
0145 detector.fromDetL.clear();
0146 detector.fromLevels.clear();
0147 fineDetectors_.emplace_back(detector);
0148 }
0149 }
0150
0151 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing: with " << detectors_.size() << " calorimetric volumes";
0152 for (unsigned int i = 0; i < detectors_.size(); i++) {
0153 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing: Calorimeter volume " << i << " " << detectors_[i].name << " LV "
0154 << detectors_[i].lv << " at level " << detectors_[i].level << " with "
0155 << detectors_[i].fromDets.size() << " neighbours";
0156 for (unsigned int k = 0; k < detectors_[i].fromDets.size(); k++)
0157 edm::LogVerbatim("CaloSim") << " Element " << k << " " << detectors_[i].fromDets[k] << " LV "
0158 << detectors_[i].fromDetL[k] << " at level " << detectors_[i].fromLevels[k];
0159 }
0160
0161 doFineCalo_ = doFineCalo_ && !(fineDetectors_.empty());
0162 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing: with " << fineDetectors_.size() << " special calorimetric volumes";
0163 for (unsigned int i = 0; i < detectors_.size(); i++)
0164 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing: Calorimeter volume " << i << " " << detectors_[i].name << " LV "
0165 << detectors_[i].lv << " at level " << detectors_[i].level;
0166 }
0167
0168 CaloTrkProcessing::~CaloTrkProcessing() {}
0169
0170 void CaloTrkProcessing::update(const BeginOfEvent* evt) { lastTrackID_ = -1; }
0171
0172 void CaloTrkProcessing::update(const G4Step* aStep) {
0173
0174
0175 G4Track* theTrack = aStep->GetTrack();
0176 int id = theTrack->GetTrackID();
0177
0178 TrackInformation* trkInfo = dynamic_cast<TrackInformation*>(theTrack->GetUserInformation());
0179
0180 if (trkInfo == nullptr) {
0181 edm::LogError("CaloSim") << "CaloTrkProcessing: No trk info !!!! abort ";
0182 throw cms::Exception("Unknown", "CaloTrkProcessing") << "cannot get trkInfo for Track " << id << "\n";
0183 }
0184
0185 if (doFineCalo_) {
0186 int prestepLV = isItCalo(aStep->GetPreStepPoint()->GetTouchable(), fineDetectors_);
0187 int poststepLV = isItCalo(aStep->GetPostStepPoint()->GetTouchable(), fineDetectors_);
0188
0189
0190 if (!trkInfo->startedInFineVolumeIsSet())
0191 trkInfo->setStartedInFineVolume(prestepLV >= 0);
0192
0193
0194 if (prestepLV < 0 && poststepLV >= 0) {
0195 #ifdef EDM_ML_DEBUG
0196 edm::LogVerbatim("DoFineCalo") << "Track " << id << " entered a fine volume:"
0197 << " pdgid=" << theTrack->GetDefinition()->GetPDGEncoding()
0198 << " theTrack->GetCurrentStepNumber()=" << theTrack->GetCurrentStepNumber()
0199 << " prestepLV=" << prestepLV << " poststepLV=" << poststepLV
0200 << " GetKineticEnergy[GeV]=" << theTrack->GetKineticEnergy() / CLHEP::GeV
0201 << " prestepPosition[cm]=("
0202 << theTrack->GetStep()->GetPreStepPoint()->GetPosition().x() / CLHEP::cm << ","
0203 << theTrack->GetStep()->GetPreStepPoint()->GetPosition().y() / CLHEP::cm << ","
0204 << theTrack->GetStep()->GetPreStepPoint()->GetPosition().z() / CLHEP::cm << ")"
0205 << " poststepPosition[cm]=("
0206 << theTrack->GetStep()->GetPostStepPoint()->GetPosition().x() / CLHEP::cm << ","
0207 << theTrack->GetStep()->GetPostStepPoint()->GetPosition().y() / CLHEP::cm << ","
0208 << theTrack->GetStep()->GetPostStepPoint()->GetPosition().z() / CLHEP::cm << ")"
0209 << " position[cm]=(" << theTrack->GetPosition().x() / CLHEP::cm << ","
0210 << theTrack->GetPosition().y() / CLHEP::cm << ","
0211 << theTrack->GetPosition().z() / CLHEP::cm << ")"
0212 << " vertex_position[cm]=(" << theTrack->GetVertexPosition().x() / CLHEP::cm << ","
0213 << theTrack->GetVertexPosition().y() / CLHEP::cm << ","
0214 << theTrack->GetVertexPosition().z() / CLHEP::cm << ")"
0215 << " GetVertexKineticEnergy[GeV]="
0216 << theTrack->GetVertexKineticEnergy() / CLHEP::GeV;
0217 #endif
0218 if (!trkInfo->startedInFineVolume() && !trkInfo->crossedBoundary()) {
0219 trkInfo->setCrossedBoundary(theTrack);
0220 #ifdef EDM_ML_DEBUG
0221 edm::LogVerbatim("DoFineCalo") << "Track " << id << " marked as boundary-crossing; sanity check:"
0222 << " theTrack->GetTrackID()=" << theTrack->GetTrackID()
0223 << " trkInfo->crossedBoundary()=" << trkInfo->crossedBoundary();
0224 #endif
0225 }
0226 #ifdef EDM_ML_DEBUG
0227 else {
0228 edm::LogVerbatim("DoFineCalo") << "Track " << id << " REENTERED a fine volume;"
0229 << " not counting this boundary crossing!";
0230 }
0231 #endif
0232
0233 }
0234 #ifdef EDM_ML_DEBUG
0235 else if (prestepLV >= 0 && poststepLV < 0) {
0236 edm::LogVerbatim("DoFineCalo") << "Track " << id << " exited a fine volume:"
0237 << " theTrack->GetCurrentStepNumber()=" << theTrack->GetCurrentStepNumber()
0238 << " prestepLV=" << prestepLV << " poststepLV=" << poststepLV
0239 << " GetKineticEnergy[GeV]=" << theTrack->GetKineticEnergy() / CLHEP::GeV
0240 << " prestepPosition[cm]=("
0241 << theTrack->GetStep()->GetPreStepPoint()->GetPosition().x() / CLHEP::cm << ","
0242 << theTrack->GetStep()->GetPreStepPoint()->GetPosition().y() / CLHEP::cm << ","
0243 << theTrack->GetStep()->GetPreStepPoint()->GetPosition().z() / CLHEP::cm << ")"
0244 << " poststepPosition[cm]=("
0245 << theTrack->GetStep()->GetPostStepPoint()->GetPosition().x() / CLHEP::cm << ","
0246 << theTrack->GetStep()->GetPostStepPoint()->GetPosition().y() / CLHEP::cm << ","
0247 << theTrack->GetStep()->GetPostStepPoint()->GetPosition().z() / CLHEP::cm << ")";
0248 }
0249 #endif
0250 }
0251
0252 if (testBeam_) {
0253 if (trkInfo->getIDonCaloSurface() == 0) {
0254 #ifdef EDM_ML_DEBUG
0255 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing set IDonCaloSurface to " << id << " at step Number "
0256 << theTrack->GetCurrentStepNumber();
0257 #endif
0258 trkInfo->setIDonCaloSurface(id, 0, 0, theTrack->GetDefinition()->GetPDGEncoding(), theTrack->GetMomentum().mag());
0259 lastTrackID_ = id;
0260 if (theTrack->GetKineticEnergy() / CLHEP::MeV > eMin_)
0261 trkInfo->putInHistory();
0262 }
0263 } else {
0264 if (putHistory_) {
0265 trkInfo->putInHistory();
0266
0267 }
0268 #ifdef EDM_ML_DEBUG
0269 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing Entered for " << id << " at stepNumber "
0270 << theTrack->GetCurrentStepNumber() << " IDonCaloSur.. "
0271 << trkInfo->getIDonCaloSurface() << " CaloCheck " << trkInfo->caloIDChecked();
0272 #endif
0273 if (trkInfo->getIDonCaloSurface() != 0) {
0274 if (trkInfo->caloIDChecked() == false) {
0275 G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
0276 const G4VTouchable* post_touch = postStepPoint->GetTouchable();
0277
0278 if (isItInside(post_touch, trkInfo->getIDCaloVolume(), trkInfo->getIDLastVolume()) > 0) {
0279 trkInfo->setIDonCaloSurface(0, -1, -1, 0, 0);
0280 } else {
0281 trkInfo->setCaloIDChecked(true);
0282 }
0283 }
0284 } else {
0285 G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
0286 const G4VTouchable* post_touch = postStepPoint->GetTouchable();
0287 int ical = isItCalo(post_touch, detectors_);
0288 if (ical >= 0) {
0289 G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0290 const G4VTouchable* pre_touch = preStepPoint->GetTouchable();
0291 int inside = isItInside(pre_touch, ical, -1);
0292 if (inside >= 0 || (theTrack->GetCurrentStepNumber() == 1)) {
0293 trkInfo->setIDonCaloSurface(
0294 id, ical, inside, theTrack->GetDefinition()->GetPDGEncoding(), theTrack->GetMomentum().mag());
0295 trkInfo->setCaloIDChecked(true);
0296 if (!doFineCalo_)
0297 trkInfo->setCrossedBoundary(theTrack);
0298 lastTrackID_ = id;
0299 if (theTrack->GetKineticEnergy() / CLHEP::MeV > eMin_)
0300 trkInfo->putInHistory();
0301 #ifdef EDM_ML_DEBUG
0302 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing: set ID on Calo " << ical << " surface (Inside " << inside
0303 << ") to " << id << " of a Track with Kinetic Energy "
0304 << theTrack->GetKineticEnergy() / CLHEP::MeV << " MeV";
0305 #endif
0306 }
0307 }
0308 }
0309 }
0310 }
0311
0312 int CaloTrkProcessing::isItCalo(const G4VTouchable* touch, const std::vector<Detector>& detectors) {
0313 int lastLevel = -1;
0314 G4LogicalVolume* lv = nullptr;
0315 for (unsigned int it = 0; it < detectors.size(); it++) {
0316 if (lastLevel != detectors[it].level) {
0317 lastLevel = detectors[it].level;
0318 lv = detLV(touch, lastLevel);
0319 #ifdef EDM_ML_DEBUG
0320 std::string name1 = "Unknown";
0321 if (lv != 0)
0322 name1 = lv->GetName();
0323 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing: volume " << name1 << " at Level " << lastLevel;
0324 int levels = detLevels(touch);
0325 if (levels > 0) {
0326 G4String name2[20];
0327 int copyno2[20];
0328 detectorLevel(touch, levels, copyno2, name2);
0329 for (int i2 = 0; i2 < levels; i2++)
0330 edm::LogVerbatim("CaloSim") << " " << i2 << " " << name2[i2] << " " << copyno2[i2];
0331 }
0332 #endif
0333 }
0334 bool ok = (lv == detectors[it].lv);
0335 if (ok)
0336 return it;
0337 }
0338 return -1;
0339 }
0340
0341 int CaloTrkProcessing::isItInside(const G4VTouchable* touch, int idcal, int idin) {
0342 int lastLevel = -1;
0343 G4LogicalVolume* lv = nullptr;
0344 int id1, id2;
0345 if (idcal < 0) {
0346 id1 = 0;
0347 id2 = static_cast<int>(detectors_.size());
0348 } else {
0349 id1 = idcal;
0350 id2 = id1 + 1;
0351 }
0352 for (int it1 = id1; it1 < id2; it1++) {
0353 if (idin < 0) {
0354 for (unsigned int it2 = 0; it2 < detectors_[it1].fromDets.size(); it2++) {
0355 if (lastLevel != detectors_[it1].fromLevels[it2]) {
0356 lastLevel = detectors_[it1].fromLevels[it2];
0357 lv = detLV(touch, lastLevel);
0358 #ifdef EDM_ML_DEBUG
0359 std::string name1 = "Unknown";
0360 if (lv != 0)
0361 name1 = lv->GetName();
0362 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing: volume " << name1 << " at Level " << lastLevel;
0363 int levels = detLevels(touch);
0364 if (levels > 0) {
0365 G4String name2[20];
0366 int copyno2[20];
0367 detectorLevel(touch, levels, copyno2, name2);
0368 for (int i2 = 0; i2 < levels; i2++)
0369 edm::LogVerbatim("CaloSim") << " " << i2 << " " << name2[i2] << " " << copyno2[i2];
0370 }
0371 #endif
0372 }
0373 bool ok = (lv == detectors_[it1].fromDetL[it2]);
0374 if (ok)
0375 return it2;
0376 }
0377 } else {
0378 lastLevel = detectors_[it1].fromLevels[idin];
0379 lv = detLV(touch, lastLevel);
0380 #ifdef EDM_ML_DEBUG
0381 std::string name1 = "Unknown";
0382 if (lv != 0)
0383 name1 = lv->GetName();
0384 edm::LogVerbatim("CaloSim") << "CaloTrkProcessing: volume " << name1 << " at Level " << lastLevel;
0385 int levels = detLevels(touch);
0386 if (levels > 0) {
0387 G4String name2[20];
0388 int copyno2[20];
0389 detectorLevel(touch, levels, copyno2, name2);
0390 for (int i2 = 0; i2 < levels; i2++)
0391 edm::LogVerbatim("CaloSim") << " " << i2 << " " << name2[i2] << " " << copyno2[i2];
0392 }
0393 #endif
0394 bool ok = (lv == detectors_[it1].fromDetL[idin]);
0395 if (ok)
0396 return idin;
0397 }
0398 }
0399 return -1;
0400 }
0401
0402 int CaloTrkProcessing::detLevels(const G4VTouchable* touch) const {
0403
0404 if (touch)
0405 return ((touch->GetHistoryDepth()) + 1);
0406 else
0407 return 0;
0408 }
0409
0410 G4LogicalVolume* CaloTrkProcessing::detLV(const G4VTouchable* touch, int currentlevel) const {
0411 G4LogicalVolume* lv = nullptr;
0412 if (touch) {
0413 int level = ((touch->GetHistoryDepth()) + 1);
0414 if (level > 0 && level >= currentlevel) {
0415 int ii = level - currentlevel;
0416 lv = touch->GetVolume(ii)->GetLogicalVolume();
0417 return lv;
0418 }
0419 }
0420 return lv;
0421 }
0422
0423 void CaloTrkProcessing::detectorLevel(const G4VTouchable* touch, int& level, int* copyno, G4String* name) const {
0424 static const std::string unknown("Unknown");
0425
0426 if (level > 0) {
0427 for (int ii = 0; ii < level; ii++) {
0428 int i = level - ii - 1;
0429 G4VPhysicalVolume* pv = touch->GetVolume(i);
0430 if (pv != nullptr)
0431 name[ii] = pv->GetName();
0432 else
0433 name[ii] = unknown;
0434 copyno[ii] = touch->GetReplicaNumber(i);
0435 }
0436 }
0437 #ifdef EDM_ML_DEBUG
0438 edm::LogVerbatim("CaloSimX") << "CaloTrkProcessing::detectorLevel "
0439 << " with " << level << ":" << detLevels(touch) << " levels";
0440 for (int ii = 0; ii < level; ii++)
0441 edm::LogVerbatim("CaloSimX") << "[" << ii << "] " << name[ii] << ":" << copyno[ii];
0442 #endif
0443 }
