File indexing completed on 2024-05-10 02:21:19
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018
0019 #include "SimG4Core/Notification/interface/TrackInformation.h"
0020 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0021 #include "SimG4Core/Physics/interface/G4ProcessTypeEnumerator.h"
0022
0023 #include "SimDataFormats/TrackingHit/interface/UpdatablePSimHit.h"
0024 #include "SimDataFormats/SimHitMaker/interface/TrackingSlaveSD.h"
0025
0026 #include "SimG4CMS/Forward/interface/TotemSD.h"
0027 #include "SimG4CMS/Forward/interface/TotemNumberMerger.h"
0028 #include "SimG4CMS/Forward/interface/TotemT1NumberingScheme.h"
0029 #include "SimG4CMS/Forward/interface/TotemT2NumberingSchemeGem.h"
0030 #include "SimG4CMS/Forward/interface/TotemRPNumberingScheme.h"
0031
0032 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0033 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0034
0035 #include "G4SDManager.hh"
0036 #include "G4Step.hh"
0037 #include "G4Track.hh"
0038 #include "G4VProcess.hh"
0039
0040 #include "G4PhysicalConstants.hh"
0041 #include <CLHEP/Units/SystemOfUnits.h>
0042
0043
0044
0045
0046
0047 TotemSD::TotemSD(const std::string& name,
0048 const SensitiveDetectorCatalog& clg,
0049 edm::ParameterSet const& p,
0050 const SimTrackManager* manager)
0051 : SensitiveTkDetector(name, clg),
0052 numberingScheme(nullptr),
0053 hcID(-1),
0054 theHC(nullptr),
0055 theManager(manager),
0056 currentHit(nullptr),
0057 theTrack(nullptr),
0058 currentPV(nullptr),
0059 unitID(0),
0060 previousUnitID(0),
0061 preStepPoint(nullptr),
0062 postStepPoint(nullptr) {
0063
0064 edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("TotemSD");
0065 int verbn = m_p.getUntrackedParameter<int>("Verbosity");
0066
0067 SetVerboseLevel(verbn);
0068
0069 slave = new TrackingSlaveSD(name);
0070
0071 if (name == "TotemHitsT1") {
0072 numberingScheme = dynamic_cast<TotemVDetectorOrganization*>(new TotemT1NumberingScheme(1));
0073 } else if (name == "TotemHitsT2Si") {
0074 numberingScheme = dynamic_cast<TotemVDetectorOrganization*>(new TotemT2NumberingSchemeGem(3));
0075 } else if (name == "TotemHitsT2Gem") {
0076 numberingScheme = dynamic_cast<TotemVDetectorOrganization*>(new TotemT2NumberingSchemeGem(4));
0077 } else if (name == "TotemHitsRP") {
0078 numberingScheme = dynamic_cast<TotemVDetectorOrganization*>(new TotemRPNumberingScheme(3));
0079 } else {
0080 edm::LogWarning("ForwardSim") << "TotemSD: ReadoutName not supported\n";
0081 }
0082 }
0083
0084 TotemSD::~TotemSD() {
0085 delete slave;
0086 delete numberingScheme;
0087 }
0088
0089 bool TotemSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
0090 getStepInfo(aStep);
0091 if (!hitExists() && edeposit > 0.) {
0092 createNewHit();
0093 return true;
0094 }
0095 if (!hitExists() && (((unitID == 1111 || unitID == 2222) && ParentId == 0 && ParticleType == 2212))) {
0096 createNewHitEvo();
0097 }
0098 return true;
0099 }
0100
0101 uint32_t TotemSD::setDetUnitId(const G4Step* aStep) {
0102 return (numberingScheme == nullptr ? 0 : numberingScheme->getUnitID(aStep));
0103 }
0104
0105 void TotemSD::Initialize(G4HCofThisEvent* HCE) {
0106 #ifdef EDM_ML_DEBUG
0107 edm::LogVerbatim("ForwardSim") << "TotemSD : Initialize called for " << GetName();
0108 #endif
0109
0110 theHC = new TotemG4HitCollection(GetName(), collectionName[0]);
0111 if (hcID < 0)
0112 hcID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
0113 HCE->AddHitsCollection(hcID, theHC);
0114
0115 tsID = -2;
0116 primID = -2;
0117 }
0118
0119 void TotemSD::EndOfEvent(G4HCofThisEvent*) {
0120
0121 int thehc_entries = theHC->entries();
0122 for (int j = 0; j < thehc_entries && j < 15000; j++) {
0123 TotemG4Hit* aHit = (*theHC)[j];
0124 #ifdef EDM_ML_DEBUG
0125 edm::LogVerbatim("ForwardSim") << "HIT NUMERO " << j << "unit ID = " << aHit->getUnitID()
0126 << "\n enrty z " << aHit->getEntry().z() << "\n theta "
0127 << aHit->getThetaAtEntry() << "\n";
0128 #endif
0129 Local3DPoint theExitPoint(0, 0, 0);
0130 Local3DPoint Entrata(aHit->getEntry().x(), aHit->getEntry().y(), aHit->getEntry().z());
0131 slave->processHits(PSimHit(Entrata,
0132 theExitPoint,
0133 aHit->getPabs(),
0134 aHit->getTof(),
0135 aHit->getEnergyLoss(),
0136 aHit->getParticleType(),
0137 aHit->getUnitID(),
0138 aHit->getTrackID(),
0139 aHit->getThetaAtEntry(),
0140 aHit->getPhiAtEntry()));
0141 }
0142 }
0143
0144 void TotemSD::PrintAll() {
0145 #ifdef EDM_ML_DEBUG
0146 edm::LogVerbatim("ForwardSim") << "TotemSD: Collection " << theHC->GetName();
0147 #endif
0148 theHC->PrintAllHits();
0149 }
0150
0151 void TotemSD::fillHits(edm::PSimHitContainer& cc, const std::string& hname) {
0152 if (slave->name() == hname) {
0153 cc = slave->hits();
0154 }
0155 }
0156
0157 void TotemSD::update(const BeginOfEvent* i) {
0158 #ifdef EDM_ML_DEBUG
0159 edm::LogVerbatim("ForwardSim") << " Dispatched BeginOfEvent for " << GetName() << " !";
0160 #endif
0161 clearHits();
0162 }
0163
0164 void TotemSD::clearHits() { slave->Initialize(); }
0165
0166 G4ThreeVector TotemSD::setToLocal(const G4ThreeVector& global) {
0167 G4ThreeVector localPoint;
0168 const G4VTouchable* touch = preStepPoint->GetTouchable();
0169 localPoint = touch->GetHistory()->GetTopTransform().TransformPoint(global);
0170 return localPoint;
0171 }
0172
0173 void TotemSD::getStepInfo(const G4Step* aStep) {
0174 preStepPoint = aStep->GetPreStepPoint();
0175 postStepPoint = aStep->GetPostStepPoint();
0176 theTrack = aStep->GetTrack();
0177 hitPoint = preStepPoint->GetPosition();
0178 currentPV = preStepPoint->GetPhysicalVolume();
0179
0180
0181 G4String name = currentPV->GetName();
0182 name.assign(name, 0, 4);
0183 G4String particleType = theTrack->GetDefinition()->GetParticleName();
0184 edeposit = aStep->GetTotalEnergyDeposit();
0185
0186 tSlice = (postStepPoint->GetGlobalTime()) / CLHEP::nanosecond;
0187 tSliceID = (int)tSlice;
0188 unitID = setDetUnitId(aStep);
0189 #ifdef EDM_ML_DEBUG
0190 edm::LogVerbatim("ForwardSim") << "UNITa " << unitID;
0191 #endif
0192 primaryID = theTrack->GetTrackID();
0193
0194 Posizio = hitPoint;
0195 Pabs = aStep->GetPreStepPoint()->GetMomentum().mag() / CLHEP::GeV;
0196 Tof = aStep->GetPostStepPoint()->GetGlobalTime() / CLHEP::nanosecond;
0197
0198 Eloss = aStep->GetTotalEnergyDeposit() / CLHEP::GeV;
0199 ParticleType = theTrack->GetDefinition()->GetPDGEncoding();
0200
0201 ThetaAtEntry = aStep->GetPreStepPoint()->GetPosition().theta() / CLHEP::deg;
0202 PhiAtEntry = aStep->GetPreStepPoint()->GetPosition().phi() / CLHEP::deg;
0203
0204 ParentId = theTrack->GetParentID();
0205 Vx = theTrack->GetVertexPosition().x();
0206 Vy = theTrack->GetVertexPosition().y();
0207 Vz = theTrack->GetVertexPosition().z();
0208 }
0209
0210 bool TotemSD::hitExists() {
0211 if (primaryID < 1) {
0212 edm::LogWarning("ForwardSim") << "***** TotemSD error: primaryID = " << primaryID << " maybe detector name changed";
0213 }
0214
0215
0216
0217 if (tSliceID == tsID && unitID == previousUnitID) {
0218 updateHit();
0219 return true;
0220 }
0221
0222
0223 if (primaryID != primID)
0224 resetForNewPrimary();
0225
0226
0227
0228
0229 bool found = false;
0230 int thehc_entries = theHC->entries();
0231 for (int j = 0; j < thehc_entries && !found; j++) {
0232 TotemG4Hit* aPreviousHit = (*theHC)[j];
0233 if (aPreviousHit->getTrackID() == primaryID && aPreviousHit->getTimeSliceID() == tSliceID &&
0234 aPreviousHit->getUnitID() == unitID) {
0235 currentHit = aPreviousHit;
0236 found = true;
0237 break;
0238 }
0239 }
0240
0241 if (found) {
0242 updateHit();
0243 return true;
0244 } else {
0245 return false;
0246 }
0247 }
0248
0249 void TotemSD::createNewHit() {
0250 #ifdef EDM_ML_DEBUG
0251 edm::LogVerbatim("ForwardSim") << "TotemSD CreateNewHit for PV " << currentPV->GetName()
0252 << " PVid = " << currentPV->GetCopyNo() << " Unit " << unitID << "\n primary "
0253 << primaryID << " time slice " << tSliceID << " For Track " << theTrack->GetTrackID()
0254 << " which is a " << theTrack->GetDefinition()->GetParticleName();
0255
0256 if (theTrack->GetTrackID() == 1) {
0257 edm::LogVerbatim("ForwardSim") << " of energy " << theTrack->GetTotalEnergy();
0258 } else {
0259 edm::LogVerbatim("ForwardSim") << " daughter of part. " << theTrack->GetParentID();
0260 }
0261
0262 if (theTrack->GetCreatorProcess() != nullptr)
0263 edm::LogVerbatim("ForwardSim") << theTrack->GetCreatorProcess()->GetProcessName();
0264 else
0265 edm::LogVerbatim("ForwardSim") << "NO process";
0266 #endif
0267
0268 currentHit = new TotemG4Hit;
0269 currentHit->setTrackID(primaryID);
0270 currentHit->setTimeSlice(tSlice);
0271 currentHit->setUnitID(unitID);
0272 currentHit->setIncidentEnergy(incidentEnergy);
0273
0274 currentHit->setPabs(Pabs);
0275 currentHit->setTof(Tof);
0276 currentHit->setEnergyLoss(Eloss);
0277 currentHit->setParticleType(ParticleType);
0278 currentHit->setThetaAtEntry(ThetaAtEntry);
0279 currentHit->setPhiAtEntry(PhiAtEntry);
0280
0281 currentHit->setEntry(Posizio.x(), Posizio.y(), Posizio.z());
0282
0283 currentHit->setParentId(ParentId);
0284 currentHit->setVx(Vx);
0285 currentHit->setVy(Vy);
0286 currentHit->setVz(Vz);
0287
0288 updateHit();
0289
0290 storeHit(currentHit);
0291 }
0292
0293 void TotemSD::createNewHitEvo() {
0294
0295
0296 currentHit = new TotemG4Hit;
0297 currentHit->setTrackID(primaryID);
0298 currentHit->setTimeSlice(tSlice);
0299 currentHit->setUnitID(unitID);
0300 currentHit->setIncidentEnergy(incidentEnergy);
0301
0302 currentHit->setPabs(Pabs);
0303 currentHit->setTof(Tof);
0304 currentHit->setEnergyLoss(Eloss);
0305 currentHit->setParticleType(ParticleType);
0306 currentHit->setThetaAtEntry(ThetaAtEntry);
0307 currentHit->setPhiAtEntry(PhiAtEntry);
0308
0309
0310
0311 currentHit->setParentId(ParentId);
0312 currentHit->setVx(Vx);
0313 currentHit->setVy(Vy);
0314 currentHit->setVz(Vz);
0315
0316 G4ThreeVector _PosizioEvo;
0317 int flagAcc = 0;
0318 _PosizioEvo = posizioEvo(Posizio, Vx, Vy, Vz, Pabs, flagAcc);
0319
0320 if (flagAcc == 1) {
0321 currentHit->setEntry(_PosizioEvo.x(), _PosizioEvo.y(), _PosizioEvo.z());
0322
0323 updateHit();
0324
0325 storeHit(currentHit);
0326 }
0327
0328 }
0329
0330 G4ThreeVector TotemSD::posizioEvo(
0331 const G4ThreeVector& Pos, double vx, double vy, double vz, double pabs, int& accettanza) {
0332 accettanza = 0;
0333
0334 G4ThreeVector PosEvo;
0335 double ThetaX = std::atan((Pos.x() - vx) / (Pos.z() - vz));
0336 double ThetaY = std::atan((Pos.y() - vy) / (Pos.z() - vz));
0337 double X_at_0 = (vx - ((Pos.x() - vx) / (Pos.z() - vz)) * vz) / 1000.;
0338 double Y_at_0 = (vy - ((Pos.y() - vy) / (Pos.z() - vz)) * vz) / 1000.;
0339
0340
0341
0342
0343
0344
0345
0346 double csi = std::abs((7000. - pabs) / 7000.);
0347
0348
0349 const int no_rp = 4;
0350 double x_par[no_rp + 1];
0351 double y_par[no_rp + 1];
0352
0353 double rp[no_rp] = {141., 149., 198., 220.};
0354
0355 double leffx[][2] = {{122.5429, -46.9312}, {125.4194, -49.1849}, {152.6, -81.157}, {98.8914, -131.8390}};
0356
0357 double leffy[][2] = {{124.2314, -55.4852}, {127.7825, -57.4503}, {179.455, -76.274}, {273.0931, -40.4626}};
0358
0359 double avx[][2] = {{0.515483, -1.0123}, {0.494122, -1.0534}, {0.2217, -1.483}, {0.004633, -1.0719}};
0360
0361 double avy[][2] = {{0.371418, -1.6327}, {0.349035, -1.6955}, {0.0815, -2.59}, {0.007592, -4.0841}};
0362
0363 double ddx[][4] = {{-0.082336, -0.092513, 112.3436, -82.5029},
0364 {-0.086927, -0.097670, 114.9513, -82.9835},
0365 {-0.092117, -0.0915, 180.6236, -82.443},
0366 {-0.050470, 0.058837, 208.1106, 20.8198}};
0367
0368 double detlim[][2] = {{0, 0}, {0.0028, 0.0021}, {0, 0}, {0.0008, 0.0013}};
0369
0370 double pipelim[][2] = {{0.026, 0.026}, {0.04, 0.04}, {0.0226, 0.0177}, {0.04, 0.04}};
0371
0372 for (int j = 0; j < no_rp; j++) {
0373
0374
0375 y_par[j] = ThetaY * (leffy[j][0] + leffy[j][1] * csi) + (avy[j][0] + avy[j][1] * csi) * Y_at_0;
0376 x_par[j] = ThetaX * (leffx[j][0] + leffx[j][1] * csi) + (avx[j][0] + avx[j][1] * csi) * X_at_0 -
0377 csi * (ddx[j][0] + (ddx[j][1] + ddx[j][2] * ThetaX) * csi + ddx[j][3] * ThetaX);
0378 }
0379
0380
0381 if (std::abs(y_par[0]) < pipelim[0][1] && sqrt((y_par[0] * y_par[0]) + (x_par[0] * x_par[0])) < pipelim[0][0]) {
0382
0383 if ((sqrt((y_par[1] * y_par[1]) + (x_par[1] * x_par[1])) < pipelim[1][0]) &&
0384 (std::abs(y_par[1]) > detlim[1][1] || x_par[1] > detlim[1][0])) {
0385 accettanza = 1;
0386 }
0387 }
0388
0389
0390 if (std::abs(y_par[0]) < pipelim[0][1] && sqrt((y_par[0]) * (y_par[0]) + (x_par[0]) * (x_par[0])) < pipelim[0][0]) {
0391
0392 if (std::abs(y_par[2]) < pipelim[2][1] && sqrt((y_par[2] * y_par[2]) + (x_par[2] * x_par[2])) < pipelim[2][0]) {
0393
0394 if ((sqrt((y_par[3] * y_par[3]) + (x_par[3] * x_par[3])) < pipelim[3][0]) &&
0395 (std::abs(y_par[3]) > detlim[3][1] || x_par[3] > detlim[3][0])) {
0396 accettanza = 1;
0397
0398 PosEvo.setX(1000 * x_par[3]);
0399 PosEvo.setY(1000 * y_par[3]);
0400 PosEvo.setZ(1000 * rp[3]);
0401 if (Pos.z() < vz)
0402 PosEvo.setZ(-1000 * rp[3]);
0403 }
0404 }
0405 }
0406
0407
0408
0409
0410
0411
0412
0413
0414
0415
0416
0417
0418
0419
0420
0421 return PosEvo;
0422 }
0423
0424 void TotemSD::updateHit() {
0425
0426 if (Eloss > 0.) {
0427 #ifdef EDM_ML_DEBUG
0428 edm::LogVerbatim("ForwardSim") << "G4TotemT1SD updateHit: add eloss " << Eloss << "\nCurrentHit=" << currentHit
0429 << ", PostStepPoint=" << postStepPoint->GetPosition();
0430 #endif
0431
0432 currentHit->setEnergyLoss(Eloss);
0433 }
0434
0435
0436 tsID = tSliceID;
0437 primID = primaryID;
0438 previousUnitID = unitID;
0439 }
0440
0441 void TotemSD::storeHit(TotemG4Hit* hit) {
0442 if (primID < 0)
0443 return;
0444 if (hit == nullptr) {
0445 edm::LogWarning("ForwardSim") << "TotemSD: hit to be stored is NULL !!";
0446 return;
0447 }
0448
0449 theHC->insert(hit);
0450 }
0451
0452 void TotemSD::resetForNewPrimary() {
0453 entrancePoint = setToLocal(hitPoint);
0454 incidentEnergy = preStepPoint->GetKineticEnergy();
0455 }