File indexing completed on 2024-05-10 02:21:19
0001
0002
0003
0004
0005
0006
0007
0008 #include "SimG4Core/Notification/interface/TrackInformation.h"
0009 #include "SimG4Core/Notification/interface/TrackInformationExtractor.h"
0010 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0011
0012 #include "SimG4CMS/Forward/interface/CastorSD.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014
0015 #include "G4SDManager.hh"
0016 #include "G4Step.hh"
0017 #include "G4Track.hh"
0018 #include "G4VProcess.hh"
0019
0020 #include "G4ios.hh"
0021 #include "G4Cerenkov.hh"
0022 #include "G4LogicalVolumeStore.hh"
0023
0024 #include <CLHEP/Units/SystemOfUnits.h>
0025 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0026 #include "Randomize.hh"
0027 #include "G4Poisson.hh"
0028
0029
0030
0031 CastorSD::CastorSD(const std::string& name,
0032 const SensitiveDetectorCatalog& clg,
0033 edm::ParameterSet const& p,
0034 const SimTrackManager* manager)
0035 : CaloSD(name, clg, p, manager),
0036 numberingScheme(nullptr),
0037 lvC3EF(nullptr),
0038 lvC3HF(nullptr),
0039 lvC4EF(nullptr),
0040 lvC4HF(nullptr),
0041 lvCAST(nullptr) {
0042 edm::ParameterSet m_CastorSD = p.getParameter<edm::ParameterSet>("CastorSD");
0043 useShowerLibrary = m_CastorSD.getParameter<bool>("useShowerLibrary");
0044 energyThresholdSL = m_CastorSD.getParameter<double>("minEnergyInGeVforUsingSLibrary");
0045 energyThresholdSL = energyThresholdSL * CLHEP::GeV;
0046
0047 non_compensation_factor = m_CastorSD.getParameter<double>("nonCompensationFactor");
0048
0049 if (useShowerLibrary) {
0050 showerLibrary = new CastorShowerLibrary(name, p);
0051 setParameterized(true);
0052 }
0053 setNumberingScheme(new CastorNumberingScheme());
0054
0055 edm::LogVerbatim("ForwardSim") << "********************************************************\n"
0056 << "* Constructing a CastorSD with name " << GetName() << "\n"
0057 << "* Using Castor Shower Library: " << useShowerLibrary << "\n"
0058 << "********************************************************";
0059
0060 const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0061 for (auto lv : *lvs) {
0062 if (strcmp((lv->GetName()).c_str(), "C3EF") == 0) {
0063 lvC3EF = lv;
0064 }
0065 if (strcmp((lv->GetName()).c_str(), "C3HF") == 0) {
0066 lvC3HF = lv;
0067 }
0068 if (strcmp((lv->GetName()).c_str(), "C4EF") == 0) {
0069 lvC4EF = lv;
0070 }
0071 if (strcmp((lv->GetName()).c_str(), "C4HF") == 0) {
0072 lvC4HF = lv;
0073 }
0074 if (strcmp((lv->GetName()).c_str(), "CAST") == 0) {
0075 lvCAST = lv;
0076 }
0077 if (lvC3EF != nullptr && lvC3HF != nullptr && lvC4EF != nullptr && lvC4HF != nullptr && lvCAST != nullptr) {
0078 break;
0079 }
0080 }
0081 edm::LogVerbatim("ForwardSim") << "CastorSD:: LogicalVolume pointers\n"
0082 << lvC3EF << " for C3EF; " << lvC3HF << " for C3HF; " << lvC4EF << " for C4EF; "
0083 << lvC4HF << " for C4HF; " << lvCAST << " for CAST. ";
0084 }
0085
0086
0087
0088 CastorSD::~CastorSD() { delete showerLibrary; }
0089
0090
0091
0092 double CastorSD::getEnergyDeposit(const G4Step* aStep) {
0093 double NCherPhot = 0.;
0094
0095
0096 auto const theTrack = aStep->GetTrack();
0097
0098
0099 auto const preStepPoint = aStep->GetPreStepPoint();
0100 auto const currentPV = preStepPoint->GetPhysicalVolume();
0101 auto const currentLV = currentPV->GetLogicalVolume();
0102
0103 #ifdef EDM_ML_DEBUG
0104 edm::LogVerbatim("ForwardSim") << "CastorSD::getEnergyDeposit:"
0105 << "\n CurrentStepNumber , TrackID , ParentID, Particle , VertexPosition ,"
0106 << " LogicalVolumeAtVertex , PV, Time"
0107 << "\n TRACKINFO: " << theTrack->GetCurrentStepNumber() << " , "
0108 << theTrack->GetTrackID() << " , " << theTrack->GetParentID() << " , "
0109 << theTrack->GetDefinition()->GetParticleName() << " , "
0110 << theTrack->GetVertexPosition() << " , "
0111 << theTrack->GetLogicalVolumeAtVertex()->GetName() << " , " << currentPV->GetName()
0112 << " , " << theTrack->GetGlobalTime();
0113 #endif
0114
0115
0116 const G4ThreeVector& hit_mom = preStepPoint->GetMomentumDirection();
0117 const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0118 double zint = hitPoint.z();
0119
0120
0121 G4int parCode = theTrack->GetDefinition()->GetPDGEncoding();
0122 bool isHad = G4TrackToParticleID::isStableHadronIon(theTrack);
0123
0124 double meanNCherPhot = 0;
0125 G4double charge = preStepPoint->GetCharge();
0126
0127 if (0.0 == charge) {
0128 return meanNCherPhot;
0129 }
0130
0131 G4double beta = preStepPoint->GetBeta();
0132 const double bThreshold = 0.67;
0133
0134 if (beta < bThreshold) {
0135 return meanNCherPhot;
0136 }
0137
0138
0139 TrackInformationExtractor TIextractor;
0140 TrackInformation& trkInfo = TIextractor(theTrack);
0141 if (!trkInfo.hasCastorHit()) {
0142 trkInfo.setCastorHitPID(parCode);
0143 }
0144 G4double stepl = aStep->GetStepLength() / CLHEP::cm;
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164 #ifdef EDM_ML_DEBUG
0165 edm::LogVerbatim("ForwardSim") << "CastorSD::getEnergyDeposit: for ID=" << theTrack->GetTrackID()
0166 << " LV: " << currentLV->GetName() << " isHad:" << isHad << " pdg=" << parCode
0167 << " castorPID=" << trkInfo.getCastorHitPID() << " sl= " << stepl
0168 << " Edep= " << aStep->GetTotalEnergyDeposit();
0169 #endif
0170 if (currentLV == lvC3EF || currentLV == lvC4EF || currentLV == lvC3HF || currentLV == lvC4HF) {
0171 const double nMedium = 1.4925;
0172
0173
0174
0175 const double photEnSpectrDE = 1.24;
0176
0177
0178
0179
0180
0181
0182
0183
0184 const double thFullRefl = 23.;
0185 const double thFullReflRad = thFullRefl * pi / 180.;
0186
0187
0188 const double thFibDir = 45.;
0189
0190
0191
0192
0193
0194 const double thFibDirRad = thFibDir * pi / 180.;
0195
0196
0197 double costh =
0198 hit_mom.z() / sqrt(hit_mom.x() * hit_mom.x() + hit_mom.y() * hit_mom.y() + hit_mom.z() * hit_mom.z());
0199 if (zint < 0)
0200 costh = -costh;
0201 double th = acos(std::min(std::max(costh, double(-1.)), double(1.)));
0202
0203
0204 if (th < 0.)
0205 th += twopi;
0206
0207
0208 double costhcher = 1. / (nMedium * beta);
0209 double thcher = acos(std::min(std::max(costhcher, double(-1.)), double(1.)));
0210
0211
0212 double DelFibPart = std::abs(th - thFibDirRad);
0213
0214
0215 double d = std::abs(tan(th) - tan(thFibDirRad));
0216
0217 double a = tan(thFibDirRad) + tan(std::abs(thFibDirRad - thFullReflRad));
0218 double r = tan(th) + tan(std::abs(th - thcher));
0219
0220
0221 double d_qz;
0222 #ifdef EDM_ML_DEBUG
0223 int variant(0);
0224 #endif
0225 if (DelFibPart > (thFullReflRad + thcher)) {
0226 d_qz = 0.;
0227 } else {
0228 if ((th + thcher) < (thFibDirRad + thFullReflRad) && (th - thcher) > (thFibDirRad - thFullReflRad)) {
0229 d_qz = 1.;
0230 #ifdef EDM_ML_DEBUG
0231 variant = 1;
0232 #endif
0233 } else {
0234 if ((thFibDirRad + thFullReflRad) < (th + thcher) && (thFibDirRad - thFullReflRad) > (th - thcher)) {
0235 d_qz = 0.;
0236 #ifdef EDM_ML_DEBUG
0237 variant = 2;
0238 #endif
0239 } else {
0240 #ifdef EDM_ML_DEBUG
0241 variant = 3;
0242 #endif
0243
0244
0245 double arg_arcos = 0.;
0246 double tan_arcos = 2. * a * d;
0247 if (tan_arcos != 0.)
0248 arg_arcos = (r * r - a * a - d * d) / tan_arcos;
0249 arg_arcos = std::abs(arg_arcos);
0250 double th_arcos = acos(std::min(std::max(arg_arcos, -1.), 1.));
0251 d_qz = std::abs(th_arcos / CLHEP::twopi);
0252 }
0253 }
0254 }
0255 #ifdef EDM_ML_DEBUG
0256 edm::LogVerbatim("ForwardSim") << "variant" << variant;
0257 #endif
0258
0259
0260
0261 meanNCherPhot = 370. * charge * charge * (1. - 1. / (nMedium * nMedium * beta * beta)) * photEnSpectrDE * stepl;
0262
0263 double scale = isHad ? non_compensation_factor : 1.0;
0264
0265 int poissNCherPhot = std::max(0, (int)G4Poisson(meanNCherPhot * scale));
0266
0267 const double effPMTandTransport = 0.19;
0268 const double ReflPower = 0.1;
0269 double proba = d_qz + (1 - d_qz) * ReflPower;
0270 NCherPhot = poissNCherPhot * effPMTandTransport * proba * 0.307;
0271 #ifdef EDM_ML_DEBUG
0272 edm::LogVerbatim("ForwardSim") << " Nph= " << NCherPhot << " Np= " << poissNCherPhot
0273 << " eff= " << effPMTandTransport << " pb= " << proba << " Nmean= " << meanNCherPhot
0274 << " q=" << charge << " beta=" << beta << " nMedium= " << nMedium << " sl= " << stepl
0275 << " Nde=" << photEnSpectrDE;
0276 #endif
0277 }
0278 return NCherPhot;
0279 }
0280
0281
0282
0283 uint32_t CastorSD::setDetUnitId(const G4Step* aStep) {
0284 return (numberingScheme == nullptr ? 0 : numberingScheme->getUnitID(aStep));
0285 }
0286
0287
0288
0289 void CastorSD::setNumberingScheme(CastorNumberingScheme* scheme) {
0290 if (scheme != nullptr) {
0291 edm::LogVerbatim("ForwardSim") << "CastorSD: updates numbering scheme for " << GetName();
0292 delete numberingScheme;
0293 numberingScheme = scheme;
0294 }
0295 }
0296
0297
0298
0299 uint32_t CastorSD::rotateUnitID(uint32_t unitID, const G4Track* track, const CastorShowerEvent& shower) {
0300
0301
0302
0303
0304
0305
0306
0307
0308 double trackPhi = track->GetPosition().phi();
0309 if (trackPhi < 0)
0310 trackPhi += 2 * M_PI;
0311
0312 double showerPhi = shower.getPrimPhi();
0313 if (showerPhi < 0)
0314 showerPhi += 2 * M_PI;
0315
0316
0317
0318 int trackOctSector = (int)(trackPhi / (M_PI / 4));
0319 int showerOctSector = (int)(showerPhi / (M_PI / 4));
0320
0321 uint32_t newUnitID;
0322 uint32_t sec = ((unitID >> 4) & 0xF);
0323 uint32_t complement = (unitID & 0xFFFFFF0F);
0324
0325
0326 double trackZ = track->GetPosition().z();
0327
0328 int aux;
0329 int dSec = 2 * (trackOctSector - showerOctSector);
0330 if (trackZ > 0)
0331 {
0332 int sec1 = sec - dSec;
0333
0334 if (sec1 < 0)
0335 sec1 += 16;
0336 if (sec1 > 15)
0337 sec1 -= 16;
0338 sec = (uint32_t)(sec1);
0339 } else {
0340 if (dSec < 0)
0341 sec += 16;
0342 sec += dSec;
0343 aux = (int)(sec / 16);
0344 sec -= aux * 16;
0345 }
0346 sec = sec << 4;
0347 newUnitID = complement | sec;
0348
0349 #ifdef EDM_ML_DEBUG
0350 if (newUnitID != unitID) {
0351 edm::LogVerbatim("ForwardSim") << "\n CastorSD::rotateUnitID: "
0352 << "\n unitID = " << unitID << "\n newUnitID = " << newUnitID;
0353 }
0354 #endif
0355
0356 return newUnitID;
0357 }
0358
0359
0360
0361 bool CastorSD::getFromLibrary(const G4Step* aStep) {
0362
0363
0364
0365
0366
0367
0368
0369
0370 auto const theTrack = aStep->GetTrack();
0371 auto parCode = theTrack->GetDefinition()->GetPDGEncoding();
0372
0373
0374 TrackInformationExtractor TIextractor;
0375 TrackInformation& trkInfo = TIextractor(theTrack);
0376 if (!trkInfo.hasCastorHit()) {
0377 trkInfo.setCastorHitPID(parCode);
0378 }
0379
0380 if (!useShowerLibrary) {
0381 return false;
0382 }
0383
0384
0385
0386 auto const preStepPoint = aStep->GetPreStepPoint();
0387 auto const currentPV = preStepPoint->GetPhysicalVolume();
0388 auto const currentLV = currentPV->GetLogicalVolume();
0389
0390 #ifdef EDM_ML_DEBUG
0391 edm::LogVerbatim("ForwardSim") << "CastorSD::getFromLibrary: for ID=" << theTrack->GetTrackID()
0392 << " parentID= " << theTrack->GetParentID() << " "
0393 << theTrack->GetDefinition()->GetParticleName() << " LV: " << currentLV->GetName()
0394 << " PV: " << currentPV->GetName() << "\n eta= " << theTrack->GetPosition().eta()
0395 << " phi= " << theTrack->GetPosition().phi()
0396 << " z(cm)= " << theTrack->GetPosition().z() / CLHEP::cm
0397 << " time(ns)= " << theTrack->GetGlobalTime()
0398 << " E(GeV)= " << theTrack->GetTotalEnergy() / CLHEP::GeV;
0399
0400 #endif
0401
0402 const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0403 const G4ThreeVector& hit_mom = preStepPoint->GetMomentumDirection();
0404 double zint = hitPoint.z();
0405 double pz = hit_mom.z();
0406
0407
0408 bool backward = (pz * zint < 0.) ? true : false;
0409
0410
0411 bool aboveThreshold = (theTrack->GetKineticEnergy() > energyThresholdSL) ? true : false;
0412
0413
0414 bool isEM = G4TrackToParticleID::isGammaElectronPositron(parCode);
0415 bool isHad = G4TrackToParticleID::isStableHadronIon(theTrack);
0416
0417
0418 double theta_max = M_PI - 3.1305;
0419 double R_mom = sqrt(hit_mom.x() * hit_mom.x() + hit_mom.y() * hit_mom.y());
0420 double theta = atan2(R_mom, std::abs(pz));
0421 bool angleok = (theta < theta_max) ? true : false;
0422
0423
0424 double R = sqrt(hitPoint.x() * hitPoint.x() + hitPoint.y() * hitPoint.y());
0425 bool dot = (zint < -14450. && R < 45.) ? true : false;
0426 bool inRange = (zint < -14700. || R > 193.) ? false : true;
0427
0428
0429 bool particleWithinShowerLibrary =
0430 aboveThreshold && (isEM || isHad) && (!backward) && inRange && !dot && angleok && currentLV == lvCAST;
0431
0432 #ifdef EDM_ML_DEBUG
0433 edm::LogVerbatim("ForwardSim") << "CastorSD::getFromLibrary: ID= " << theTrack->GetTrackID() << " E>E0 "
0434 << aboveThreshold << " isEM " << isEM << " isHad " << isHad << " backword " << backward
0435 << " Ok " << (inRange && !dot) << " angle " << angleok
0436 << " LV: " << currentLV->GetName() << " " << (currentLV == lvCAST) << " "
0437 << particleWithinShowerLibrary << " Edep= " << aStep->GetTotalEnergyDeposit();
0438 #endif
0439
0440
0441
0442 if (!particleWithinShowerLibrary) {
0443 return false;
0444 }
0445
0446
0447
0448 bool isKilled(true);
0449 CastorShowerEvent hits = showerLibrary->getShowerHits(aStep, isKilled);
0450
0451 int primaryID = setTrackID(aStep);
0452
0453
0454 resetForNewPrimary(aStep);
0455
0456 #ifdef EDM_ML_DEBUG
0457 edm::LogVerbatim("ForwardSim") << "CastorSD::getFromLibrary: " << hits.getNhit() << " hits for " << GetName()
0458 << " from " << theTrack->GetDefinition()->GetParticleName() << " of "
0459 << preStepPoint->GetKineticEnergy() / CLHEP::GeV << " GeV and trackID "
0460 << theTrack->GetTrackID() << " isHad: " << isHad;
0461 #endif
0462
0463
0464 double E_track = preStepPoint->GetTotalEnergy();
0465 double E_SLhit = hits.getPrimE() * CLHEP::GeV;
0466 double scale = E_track / E_SLhit;
0467
0468
0469 if (isHad) {
0470 scale *= non_compensation_factor;
0471 }
0472
0473 for (unsigned int i = 0; i < hits.getNhit(); ++i) {
0474
0475 double nPhotoElectrons = hits.getNphotons(i) * scale;
0476
0477 if (isEM) {
0478 edepositEM = nPhotoElectrons;
0479 edepositHAD = 0.;
0480 } else {
0481 edepositEM = 0.;
0482 edepositHAD = nPhotoElectrons;
0483 }
0484
0485
0486 double time = hits.getTime(i);
0487
0488
0489 unsigned int unitID = hits.getDetID(i);
0490
0491
0492
0493
0494 unsigned int rotatedUnitID = rotateUnitID(unitID, theTrack, hits);
0495 currentID[0].setID(rotatedUnitID, time, primaryID, 0);
0496 processHit(aStep);
0497 }
0498 return isKilled;
0499 }