15 #include "G4EventManager.hh"
16 #include "G4LogicalVolumeStore.hh"
17 #include "G4LogicalVolume.hh"
18 #include "G4SDManager.hh"
21 #include "G4VProcess.hh"
22 #include "G4GFlashSpot.hh"
23 #include "G4ParticleTable.hh"
25 #include "G4SystemOfUnits.hh"
26 #include "G4PhysicalConstants.hh"
41 G4VGFlashSensitiveDetector(),
44 m_trackManager(manager),
46 ignoreTrackID(ignoreTkID),
54 std::vector<double> eminHits = m_CaloSD.
getParameter<std::vector<double>>(
"EminHits");
55 std::vector<double> tmaxHits = m_CaloSD.
getParameter<std::vector<double>>(
"TmaxHits");
56 std::vector<std::string> hcn = m_CaloSD.
getParameter<std::vector<std::string>>(
"HCNames");
57 std::vector<int> useResMap = m_CaloSD.
getParameter<std::vector<int>>(
"UseResponseTables");
58 std::vector<double> eminHitX = m_CaloSD.
getParameter<std::vector<double>>(
"EminHitsDepth");
67 double beamZ = m_CaloSD.
getParameter<
double>(
"BeamPosition") * CLHEP::cm;
68 correctT = beamZ / CLHEP::c_light / CLHEP::nanosecond;
71 SetVerboseLevel(verbn);
73 for (
unsigned int k = 0;
k < hcn.size(); ++
k) {
75 if (
k < eminHits.size())
77 if (
k < eminHitX.size())
79 if (
k < tmaxHits.size())
81 if (
k < useResMap.size() && useResMap[
k] > 0) {
87 slave = std::make_unique<CaloSlaveSD>(
name);
103 <<
" before saving the hit\n Correct TOF globally by " <<
correctT
105 <<
" ns and if energy is above " <<
eminHit /
CLHEP::MeV <<
" MeV (for depth 0) or "
119 const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
120 std::vector<G4LogicalVolume*>::const_iterator lvcite;
122 G4LogicalVolume* lv =
nullptr;
124 for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
125 if ((*lvcite)->GetName() ==
name) {
143 <<
" pointer to LV: " <<
detector.lv;
147 edm::LogError(
"CaloSim") <<
"CaloSD: Cannot find CaloSimulationParameters";
148 throw cms::Exception(
"Unknown",
"CaloSD") <<
"Cannot find CaloSimulationParameters\n";
159 edm::LogVerbatim(
"CaloSim") <<
"CaloSD::" << GetName() <<
" ID= " << aStep->GetTrack()->GetTrackID()
160 <<
" prID= " << aStep->GetTrack()->GetParentID()
161 <<
" Eprestep= " << aStep->GetPreStepPoint()->GetKineticEnergy()
162 <<
" step= " << aStep->GetStepLength() <<
" Edep= " << aStep->GetTotalEnergyDeposit();
168 aStep->GetTrack()->SetTrackStatus(fStopAndKill);
169 auto tv = aStep->GetSecondary();
170 auto vol = aStep->GetPreStepPoint()->GetPhysicalVolume();
171 for (
auto& tk : *tv) {
172 if (tk->GetVolume() == vol) {
173 tk->SetTrackStatus(fStopAndKill);
183 auto const theTrack = aStep->GetTrack();
186 double time = theTrack->GetGlobalTime() / nanosecond;
191 if (aStep->GetTotalEnergyDeposit() > 0.0 && (!
ignoreReject)) {
192 const G4TouchableHistory* touch = static_cast<const G4TouchableHistory*>(theTrack->GetTouchable());
194 <<
" Detector: " << GetName() <<
" trackID= " << theTrack->GetTrackID() <<
" "
195 << theTrack->GetDefinition()->GetParticleName()
196 <<
"\n Edep= " << aStep->GetTotalEnergyDeposit()
197 <<
" PV: " << touch->GetVolume(0)->GetName()
198 <<
" PVid= " << touch->GetReplicaNumber(0) <<
" MVid= " << touch->GetReplicaNumber(1);
203 if (aStep->GetTotalEnergyDeposit() == 0.0) {
215 G4TouchableHistory* touch = (G4TouchableHistory*)(theTrack->GetTouchable());
216 edm::LogVerbatim(
"CaloSim") <<
"CaloSD::" << GetName() <<
" PV:" << touch->GetVolume(0)->GetName()
217 <<
" PVid=" << touch->GetReplicaNumber(0) <<
" MVid=" << touch->GetReplicaNumber(1)
219 <<
edepositHAD <<
" ID=" << theTrack->GetTrackID() <<
" pID=" << theTrack->GetParentID()
220 <<
" E=" << theTrack->GetKineticEnergy() <<
" S=" << aStep->GetStepLength() <<
"\n "
221 << theTrack->GetDefinition()->GetParticleName() <<
" primaryID= " << primaryID
234 const G4Track*
track = aSpot->GetOriginatorTrack()->GetPrimaryTrack();
236 double edep = aSpot->GetEnergySpot()->GetEnergy();
242 G4StepPoint* fFakePreStepPoint = fFakeStep.GetPreStepPoint();
243 G4StepPoint* fFakePostStepPoint = fFakeStep.GetPostStepPoint();
244 fFakePreStepPoint->SetPosition(aSpot->GetPosition());
245 fFakePostStepPoint->SetPosition(aSpot->GetPosition());
247 G4TouchableHandle fTouchableHandle = aSpot->GetTouchableHandle();
248 fFakePreStepPoint->SetTouchableHandle(fTouchableHandle);
249 fFakeStep.SetTotalEnergyDeposit(edep);
276 posGlobal = aSpot->GetEnergySpot()->GetPosition();
305 int level = ((touch->GetHistoryDepth()) + 1);
309 G4LogicalVolume* lv = touch->GetVolume(
ii)->GetLogicalVolume();
312 std::string name1 = (lv == 0) ?
"Unknown" : lv->GetName();
327 edm::LogVerbatim(
"CaloSim") <<
"CaloSD : Initialize called for " << GetName();
359 theHC->PrintAllHits();
364 edm::LogVerbatim(
"CaloSim") <<
"CaloSD: Tries to transfer " <<
slave.get()->hits().size() <<
" hits for "
365 <<
slave.get()->name() <<
" " << hname;
367 if (
slave.get()->name() == hname) {
370 slave.get()->Clean();
374 return touch->GetHistory()->GetTopTransform().TransformPoint(global);
378 return touch->GetHistory()->GetTopTransform().Inverse().TransformPoint(
local);
389 posGlobal = aStep->GetPreStepPoint()->GetPosition();
400 std::map<CaloHitID, CaloG4Hit*>::const_iterator it =
hitMap.find(
currentID);
408 int maxhit =
nhits - 1;
410 for (
int j = maxhit;
j > minhit; --
j) {
433 << theTrack->GetDefinition()->GetParticleName()
434 <<
" E(GeV)= " << theTrack->GetKineticEnergy() /
CLHEP::GeV
435 <<
" parentID= " << theTrack->GetParentID() <<
"\n Ein= " <<
incidentEnergy
461 etrack = theTrack->GetKineticEnergy();
477 if (trkh !=
nullptr) {
505 auto const preStepPoint = aStep->GetPreStepPoint();
510 edm::LogVerbatim(
"CaloSim") <<
"CaloSD::resetForNewPrimary for " << GetName()
519 double charge = aStep->GetPreStepPoint()->GetCharge();
520 double length = aStep->GetStepLength();
522 if (
charge != 0. && length > 0.) {
523 double density = aStep->GetPreStepPoint()->GetMaterial()->GetDensity();
524 double dedx = aStep->GetTotalEnergyDeposit() / length;
526 double c = birk2 * rkb * rkb;
529 weight = 1. / (1. + rkb * dedx +
c * dedx * dedx);
531 edm::LogVerbatim(
"CaloSim") <<
"CaloSD::getAttenuation in " << aStep->GetPreStepPoint()->GetMaterial()->GetName()
532 <<
" Charge " <<
charge <<
" dE/dx " << dedx <<
" Birk Const " << rkb <<
", " <<
c
533 <<
" Weight = " <<
weight <<
" dE " << aStep->GetTotalEnergyDeposit();
543 edm::LogVerbatim(
"CaloSim") <<
"CaloSD: Dispatched BeginOfEvent for " << GetName() <<
" !";
550 int id = (*trk)()->GetTrackID();
552 int lastTrackID = -1;
555 if (
id == lastTrackID) {
557 if (trksForThisEvent !=
nullptr) {
558 int it = (
int)(trksForThisEvent->size()) - 1;
564 edm::LogVerbatim(
"CaloSim") <<
"CaloSD: get track " << it <<
" from Container of size "
565 << trksForThisEvent->size() <<
" with ID " << trkH->
trackID();
567 edm::LogVerbatim(
"CaloSim") <<
"CaloSD: get track " << it <<
" from Container of size "
568 << trksForThisEvent->size() <<
" with no ID";
589 int hc_entries =
theHC->entries();
590 for (
int i = 0;
i < hc_entries; ++
i) {
595 double x = (*theHC)[
i]->getEM();
598 x = (*theHC)[
i]->getHadr();
601 tt += (*theHC)[
i]->getTimeSlice();
602 ee += (*theHC)[
i]->getIncidentEnergy();
607 double norm = (
count > 0) ? 1.0 /
count : 0.0;
622 <<
" hits not passing cuts\n EmeanEM= " << eEM <<
" ErmsEM= " << eEM2
623 <<
"\n EmeanHAD= " << eHAD <<
" ErmsHAD= " << eHAD2 <<
" TimeMean= " <<
tt
624 <<
" E0mean= " << ee <<
" Zglob= " << zglob <<
" Zloc= " << zloc <<
" ";
627 std::vector<std::unique_ptr<CaloG4Hit>>().
swap(
reusehit);
637 edm::LogVerbatim(
"CaloSim") <<
"CaloSD: Clears hit vector for " << GetName()
638 <<
" and initialise slave: " <<
slave.get()->name();
640 slave.get()->Initialize();
663 <<
":" << aTrack->GetTrackID() <<
":" << primaryID;
666 primaryID = aTrack->GetTrackID();
668 edm::LogWarning(
"CaloSim") <<
"CaloSD: Problem with primaryID **** set by force to TkID **** " << primaryID;
675 auto const theTrack = aStep->GetTrack();
678 if (primaryID <= 0) {
679 primaryID = theTrack->GetTrackID();
683 << theTrack->GetTrackID() <<
":" << primaryID;
691 <<
" trackID= " << aStep->GetTrack()->GetTrackID() <<
" primaryID= " << primaryID;
700 if (
hit->getDepth() > 0)
703 edm::LogVerbatim(
"CaloSim") <<
"CaloSD::filterHit(..) Depth " <<
hit->getDepth() <<
" Emin = " << emin <<
" ("
750 <<
" by SimTrackManager Status " <<
ok;
755 slave.get()->processHits(
759 << aHit->
getDepth() <<
" due to " << tkID <<
" in time " <<
time <<
" of energy "
770 primary = (*trk)()->GetTrackID();
784 if (
theHC->entries() > 0)
790 std::vector<CaloG4Hit*>* theCollection =
theHC->GetVector();
800 std::vector<CaloG4Hit*> hitvec;
802 hitvec.swap(*theCollection);
805 edm::LogVerbatim(
"CaloSim") <<
"CaloSD::cleanHitCollection: sort hits in buffer starting from "
807 for (
unsigned int i = 0;
i < hitvec.size(); ++
i)
813 for (
unsigned int j =
i + 1;
j < hitvec.size() &&
equal(hitvec[
i], hitvec[
j]); ++
j) {
816 (*hitvec[
i]).addEnergyDeposit(*hitvec[
j]);
817 (*hitvec[
j]).setEM(0.);
818 (*hitvec[
j]).setHadr(0.);
825 edm::LogVerbatim(
"CaloSim") <<
"CaloSD: cleanHitCollection merge the hits in buffer ";
826 for (
unsigned int i = 0;
i < hitvec.size(); ++
i)
831 std::stable_partition(hitvec.begin() +
cleanIndex, hitvec.end(), [](
CaloG4Hit*
p) {
return p !=
nullptr; }),
834 edm::LogVerbatim(
"CaloSim") <<
"CaloSD::cleanHitCollection: remove the merged hits in buffer,"
835 <<
" new size = " << hitvec.size();
836 for (
unsigned int i = 0;
i < hitvec.size(); ++
i)
839 hitvec.swap(*theCollection);
845 <<
" Size of reusehit= " <<
reusehit.size()
846 <<
"\n starting hit selection from index = " <<
cleanIndex;
850 for (
unsigned int i =
cleanIndex;
i < theCollection->size(); ++
i) {
866 reusehit.emplace_back((*theCollection)[
i]);
867 (*theCollection)[
i] =
nullptr;
874 <<
" Number of added hit = " << addhit;
879 for (
int ii = addhit - 1;
ii >= 0; --
ii) {
887 theCollection->erase(
888 std::stable_partition(
889 theCollection->begin() +
cleanIndex, theCollection->end(), [](
CaloG4Hit*
p) {
return p !=
nullptr; }),
890 theCollection->end());
892 edm::LogVerbatim(
"CaloSim") <<
"CaloSD: hit collection after selection, size = " <<
theHC->entries();
893 theHC->PrintAllHits();