Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-05-01 02:03:30

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