Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-09 03:15:26

0001 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0002 #include "SimG4Core/Notification/interface/TrackWithHistory.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 #include "G4SystemOfUnits.hh"
0016 #include "DD4hep/Filter.h"
0017 
0018 #include <sstream>
0019 // #define EDM_ML_DEBUG
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   //Initialise the parameter set
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   // Debug prints
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(static_cast<std::string>(dd4hep::dd::noNamespace((*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(static_cast<std::string>(dd4hep::dd::noNamespace((*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(static_cast<std::string>(dd4hep::dd::noNamespace((*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   // define if you are at the surface of CALO
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     // Once per track, determine whether track started in fine volume
0190     if (!trkInfo->startedInFineVolumeIsSet())
0191       trkInfo->setStartedInFineVolume(prestepLV >= 0);
0192 
0193     // Boundary-crossing logic
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       //      trkInfo->setAncestor();
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   //Return number of levels
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   //Get name and copy numbers
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 }