CMS 3D CMS Logo
ZdcSD.cc
Go to the documentation of this file.
1 // File: ZdcSD.cc
3 // Date: 03.01
4 // Description: Sensitive Detector class for Zdc
5 // Modifications:
7 #include <memory>
8 
9 #include "SimG4CMS/Forward/interface/ZdcSD.h"
10 #include "FWCore/MessageLogger/interface/MessageLogger.h"
11 #include "FWCore/ParameterSet/interface/ParameterSet.h"
12 #include "Geometry/Records/interface/IdealGeometryRecord.h"
13 #include "SimG4Core/Notification/interface/TrackInformation.h"
14 
15 #include "G4SDManager.hh"
16 #include "G4Step.hh"
17 #include "G4Track.hh"
18 #include "G4VProcess.hh"
19 #include "G4ios.hh"
20 #include "G4Cerenkov.hh"
21 #include "G4ParticleTable.hh"
22 #include "CLHEP/Units/GlobalSystemOfUnits.h"
23 #include "CLHEP/Units/GlobalPhysicalConstants.h"
24 #include "Randomize.hh"
25 #include "G4Poisson.hh"
26 
27 //#define EDM_ML_DEBUG
28 
29 ZdcSD::ZdcSD(const std::string& name,
30  const SensitiveDetectorCatalog& clg,
31  edm::ParameterSet const& p,
32  const SimTrackManager* manager)
33  : CaloSD(name, clg, p, manager) {
34  edm::ParameterSet m_ZdcSD = p.getParameter<edm::ParameterSet>("ZdcSD");
35  useShowerLibrary = m_ZdcSD.getParameter<bool>("UseShowerLibrary");
36  useShowerHits = m_ZdcSD.getParameter<bool>("UseShowerHits");
37  zdcHitEnergyCut = m_ZdcSD.getParameter<double>("ZdcHitEnergyCut") * GeV;
38  thFibDir = m_ZdcSD.getParameter<double>("FiberDirection");
39  verbosity = m_ZdcSD.getParameter<int>("Verbosity");
40  int verbn = verbosity / 10;
41  verbosity %= 10;
42  setNumberingScheme(new ZdcNumberingScheme(verbn));
43 
44  edm::LogVerbatim("ForwardSim") << "***************************************************\n"
45  << "* *\n"
46  << "* Constructing a ZdcSD with name " << name << " *\n"
47  << "* *\n"
48  << "***************************************************";
49 
50  edm::LogVerbatim("ForwardSim") << "\nUse of shower library is set to " << useShowerLibrary
51  << "\nUse of Shower hits method is set to " << useShowerHits;
52 
53  edm::LogVerbatim("ForwardSim") << "\nEnergy Threshold Cut set to " << zdcHitEnergyCut / GeV << " (GeV)";
54 
55  if (useShowerLibrary) {
56  showerLibrary = std::make_unique<ZdcShowerLibrary>(name, p);
57  setParameterized(true);
58  } else {
59  showerLibrary.reset(nullptr);
60  }
61 }
62 
63 void ZdcSD::initRun() { hits.clear(); }
64 
65 bool ZdcSD::getFromLibrary(const G4Step* aStep) {
66  bool ok = true;
67 
68  auto const preStepPoint = aStep->GetPreStepPoint();
69 
70  double etrack = preStepPoint->GetKineticEnergy();
71  int primaryID = setTrackID(aStep);
72 
73  hits.clear();
74 
75  // Reset entry point for new primary
76  resetForNewPrimary(aStep);
77 
78  if (etrack >= zdcHitEnergyCut) {
79  // create hits only if above threshold
80 
81 #ifdef EDM_ML_DEBUG
82  auto const theTrack = aStep->GetTrack();
83  edm::LogVerbatim("ForwardSim") << "----------------New track------------------------------\n"
84  << "Incident EnergyTrack: " << etrack << " MeV \n"
85  << "Zdc Cut Energy for Hits: " << zdcHitEnergyCut << " MeV \n"
86  << "ZdcSD::getFromLibrary " << hits.size() << " hits for " << GetName() << " of "
87  << primaryID << " with " << theTrack->GetDefinition()->GetParticleName() << " of "
88  << etrack << " MeV\n";
89 #endif
90  hits.swap(showerLibrary.get()->getHits(aStep, ok));
91  }
92 
93  incidentEnergy = etrack;
94  entrancePoint = preStepPoint->GetPosition();
95  for (unsigned int i = 0; i < hits.size(); i++) {
96  posGlobal = hits[i].position;
97  entranceLocal = hits[i].entryLocal;
98  double time = hits[i].time;
99  unsigned int unitID = hits[i].detID;
100  edepositHAD = hits[i].DeHad;
101  edepositEM = hits[i].DeEM;
102  currentID.setID(unitID, time, primaryID, 0);
103  processHit(aStep);
104 
105 #ifdef EDM_ML_DEBUG
106  edm::LogVerbatim("ForwardSim") << "ZdcSD: Final Hit number:" << i << "-->"
107  << "New HitID: " << currentHit->getUnitID()
108  << " New Hit trackID: " << currentHit->getTrackID()
109  << " New EM Energy: " << currentHit->getEM() / GeV
110  << " New HAD Energy: " << currentHit->getHadr() / GeV
111  << " New HitEntryPoint: " << currentHit->getEntryLocal()
112  << " New IncidentEnergy: " << currentHit->getIncidentEnergy() / GeV
113  << " New HitPosition: " << posGlobal;
114 #endif
115  }
116  return ok;
117 }
118 
119 double ZdcSD::getEnergyDeposit(const G4Step* aStep) {
120  double NCherPhot = 0.;
121 
122  // preStepPoint information
123  G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
124  G4VPhysicalVolume* currentPV = preStepPoint->GetPhysicalVolume();
125  const G4String& nameVolume = currentPV->GetName();
126 
127  const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
128  const G4ThreeVector& hit_mom = preStepPoint->GetMomentumDirection();
129  G4double stepL = aStep->GetStepLength() / cm;
130  G4double beta = preStepPoint->GetBeta();
131  G4double charge = preStepPoint->GetCharge();
132 
133  // theTrack information
134  G4Track* theTrack = aStep->GetTrack();
135  G4String particleType = theTrack->GetDefinition()->GetParticleName();
136  G4ThreeVector localPoint = theTrack->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
137 
138 #ifdef EDM_ML_DEBUG
139  const G4ThreeVector& vert_mom = theTrack->GetVertexMomentumDirection();
140 
141  // calculations
142  float costheta =
143  vert_mom.z() / sqrt(vert_mom.x() * vert_mom.x() + vert_mom.y() * vert_mom.y() + vert_mom.z() * vert_mom.z());
144  float theta = std::acos(std::min(std::max(costheta, -1.f), 1.f));
145  float eta = -std::log(std::tan(theta * 0.5f));
146  float phi = -100.;
147  if (vert_mom.x() != 0)
148  phi = std::atan2(vert_mom.y(), vert_mom.x());
149  if (phi < 0.)
150  phi += twopi;
151 
152  // Get the total energy deposit
153  double stepE = aStep->GetTotalEnergyDeposit();
154 
155  // postStepPoint information
156  G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
157  G4VPhysicalVolume* postPV = postStepPoint->GetPhysicalVolume();
158  const G4String& postnameVolume = postPV->GetName();
159  edm::LogVerbatim("ForwardSim") << "ZdcSD:: getEnergyDeposit: \n"
160  << " preStepPoint: " << nameVolume << "," << stepL << "," << stepE << "," << beta
161  << "," << charge << "\n"
162  << " postStepPoint: " << postnameVolume << "," << costheta << "," << theta << ","
163  << eta << "," << phi << "," << particleType << " id= " << theTrack->GetTrackID()
164  << " Etot(GeV)= " << theTrack->GetTotalEnergy() / GeV;
165 #endif
166  const double bThreshold = 0.67;
167  if ((beta > bThreshold) && (charge != 0) && (nameVolume == "ZDC_EMFiber" || nameVolume == "ZDC_HadFiber")) {
168 #ifdef EDM_ML_DEBUG
169  edm::LogVerbatim("ForwardSim") << "ZdcSD:: getEnergyDeposit: pass ";
170 #endif
171  const float nMedium = 1.4925;
172  // float photEnSpectrDL = 10714.285714;
173  // photEnSpectrDL = (1./400.nm-1./700.nm)*10000000.cm/nm; /* cm-1 */
174 
175  const float photEnSpectrDE = 1.24;
176  // E = 2pi*(1./137.)*(eV*cm/370.)/lambda = 12.389184*(eV*cm)/lambda
177  // Emax = 12.389184*(eV*cm)/400nm*10-7cm/nm = 3.01 eV
178  // Emin = 12.389184*(eV*cm)/700nm*10-7cm/nm = 1.77 eV
179  // delE = Emax - Emin = 1.24 eV
180 
181  const float effPMTandTransport = 0.15;
182 
183  // Check these values
184  const float thFullRefl = 23.;
185  float thFullReflRad = thFullRefl * pi / 180.;
186 
187  float thFibDirRad = thFibDir * pi / 180.;
188 
189  // at which theta the point is located:
190  // float th1 = hitPoint.theta();
191 
192  // theta of charged particle in LabRF(hit momentum direction):
193  float costh = hit_mom.z() / sqrt(hit_mom.x() * hit_mom.x() + hit_mom.y() * hit_mom.y() + hit_mom.z() * hit_mom.z());
194  float th = acos(std::min(std::max(costh, -1.f), 1.f));
195  // just in case (can do both standard ranges of phi):
196  if (th < 0.)
197  th += twopi;
198 
199  // theta of cone with Cherenkov photons w.r.t.direction of charged part.:
200  float costhcher = 1. / (nMedium * beta);
201  float thcher = acos(std::min(std::max(costhcher, -1.f), 1.f));
202 
203  // diff thetas of charged part. and quartz direction in LabRF:
204  float DelFibPart = std::abs(th - thFibDirRad);
205 
206  // define real distances:
207  float d = std::abs(std::tan(th) - std::tan(thFibDirRad));
208 
209  float a = std::tan(thFibDirRad) + std::tan(std::abs(thFibDirRad - thFullReflRad));
210  float r = std::tan(th) + std::tan(std::abs(th - thcher));
211 
212  // define losses d_qz in cone of full reflection inside quartz direction
213  float d_qz = -1;
214 #ifdef EDM_ML_DEBUG
215  float variant = -1;
216 #endif
217  // if (d > (r+a))
218  if (DelFibPart > (thFullReflRad + thcher)) {
219 #ifdef EDM_ML_DEBUG
220  variant = 0.;
221 #endif
222  d_qz = 0.;
223  } else {
224  // if ((DelFibPart + thcher) < thFullReflRad ) [(d+r) < a]
225  if ((th + thcher) < (thFibDirRad + thFullReflRad) && (th - thcher) > (thFibDirRad - thFullReflRad)) {
226 #ifdef EDM_ML_DEBUG
227  variant = 1.;
228 #endif
229  d_qz = 1.;
230  } else {
231  // if ((thcher - DelFibPart ) > thFullReflRad ) [(r-d) > a]
232  if ((thFibDirRad + thFullReflRad) < (th + thcher) && (thFibDirRad - thFullReflRad) > (th - thcher)) {
233 #ifdef EDM_ML_DEBUG
234  variant = 2.;
235 #endif
236  d_qz = 0.;
237  } else {
238 #ifdef EDM_ML_DEBUG
239  variant = 3.; // d_qz is calculated below
240 #endif
241  // use crossed length of circles(cone projection) - dC1/dC2 :
242  float arg_arcos = 0.;
243  float tan_arcos = 2. * a * d;
244  if (tan_arcos != 0.)
245  arg_arcos = (r * r - a * a - d * d) / tan_arcos;
246  // std::cout.testOut << " d_qz: " << r << "," << a << "," << d << " " << tan_arcos << " " << arg_arcos;
247  arg_arcos = std::abs(arg_arcos);
248  // std::cout.testOut << "," << arg_arcos;
249  float th_arcos = acos(std::min(std::max(arg_arcos, -1.f), 1.f));
250  // std::cout.testOut << " " << th_arcos;
251  d_qz = th_arcos / twopi;
252  // std::cout.testOut << " " << d_qz;
253  d_qz = std::abs(d_qz);
254  // std::cout.testOut << "," << d_qz;
255  }
256  }
257  }
258  double meanNCherPhot = 0.;
259  int poissNCherPhot = 0;
260  if (d_qz > 0) {
261  meanNCherPhot = 370. * charge * charge * (1. - 1. / (nMedium * nMedium * beta * beta)) * photEnSpectrDE * stepL;
262 
263  poissNCherPhot = std::max((int)G4Poisson(meanNCherPhot), 0);
264  NCherPhot = poissNCherPhot * effPMTandTransport * d_qz;
265  }
266 
267 #ifdef EDM_ML_DEBUG
268  edm::LogVerbatim("ForwardSim") << "ZdcSD:: getEnergyDeposit: gED: " << stepE << "," << costh << "," << th << ","
269  << costhcher << "," << thcher << "," << DelFibPart << "," << d << "," << a << ","
270  << r << "," << hitPoint << "," << hit_mom << "," << vert_mom << "," << localPoint
271  << "," << charge << "," << beta << "," << stepL << "," << d_qz << "," << variant
272  << "," << meanNCherPhot << "," << poissNCherPhot << "," << NCherPhot;
273 #endif
274  // --constants-----------------
275  // << "," << photEnSpectrDE
276  // << "," << nMedium
277  // << "," << bThreshold
278  // << "," << thFibDirRad
279  // << "," << thFullReflRad
280  // << "," << effPMTandTransport
281  // --other variables-----------
282  // << "," << curprocess
283  // << "," << nameProcess
284  // << "," << name
285  // << "," << rad
286  // << "," << mat
287 
288  } else {
289  // determine failure mode: beta, charge, and/or nameVolume
290  if (beta <= bThreshold)
291  edm::LogVerbatim("ForwardSim") << "ZdcSD:: getEnergyDeposit: fail beta=" << beta;
292  if (charge == 0)
293  edm::LogVerbatim("ForwardSim") << "ZdcSD:: getEnergyDeposit: fail charge=0";
294  if (!(nameVolume == "ZDC_EMFiber" || nameVolume == "ZDC_HadFiber"))
295  edm::LogVerbatim("ForwardSim") << "ZdcSD:: getEnergyDeposit: fail nv=" << nameVolume;
296  }
297 
298  return NCherPhot;
299 }
300 
301 uint32_t ZdcSD::setDetUnitId(const G4Step* aStep) {
302  return (numberingScheme.get() == nullptr ? 0 : numberingScheme.get()->getUnitID(aStep));
303 }
304 
305 void ZdcSD::setNumberingScheme(ZdcNumberingScheme* scheme) {
306  if (scheme != nullptr) {
307  edm::LogVerbatim("ForwardSim") << "ZdcSD: updates numbering scheme for " << GetName();
308  numberingScheme.reset(scheme);
309  }
310 }
CaloSD::edepositHAD
float edepositHAD
Definition: CaloSD.h:137
SimTrackManager
Definition: SimTrackManager.h:35
CaloG4Hit::getTrackID
int getTrackID() const
Definition: CaloG4Hit.h:64
mps_fire.i
i
Definition: mps_fire.py:428
ZdcSD::setNumberingScheme
void setNumberingScheme(ZdcNumberingScheme *scheme)
Definition: ZdcSD.cc:305
MessageLogger.h
ZdcSD::thFibDir
double thFibDir
Definition: ZdcSD.h:33
CaloSD::currentHit
CaloG4Hit * currentHit
Definition: CaloSD.h:143
ZdcSD::numberingScheme
std::unique_ptr< ZdcNumberingScheme > numberingScheme
Definition: ZdcSD.h:37
ZdcSD::hits
std::vector< ZdcShowerLibrary::Hit > hits
Definition: ZdcSD.h:38
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
ZdcSD::useShowerHits
bool useShowerHits
Definition: ZdcSD.h:32
min
T min(T a, T b)
Definition: MathUtil.h:58
HLT_FULL_cff.beta
beta
Definition: HLT_FULL_cff.py:8651
CaloG4Hit::getUnitID
uint32_t getUnitID() const
Definition: CaloG4Hit.h:66
protons_cff.time
time
Definition: protons_cff.py:35
CaloSD::edepositEM
float edepositEM
Definition: CaloSD.h:137
CaloG4Hit::getEntryLocal
math::XYZPoint getEntryLocal() const
Definition: CaloG4Hit.h:49
ZdcSD.h
convertSQLiteXML.ok
bool ok
Definition: convertSQLiteXML.py:98
ZdcSD::verbosity
int verbosity
Definition: ZdcSD.h:31
CaloSD::setTrackID
virtual int setTrackID(const G4Step *)
Definition: CaloSD.cc:839
CaloHitID::setID
void setID(uint32_t unitID, double timeSlice, int trackID, uint16_t depth=0)
Definition: CaloHitID.cc:41
PVValHelper::eta
Definition: PVValidationHelpers.h:70
ZdcNumberingScheme
Definition: ZdcNumberingScheme.h:13
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
CaloSD::entrancePoint
G4ThreeVector entrancePoint
Definition: CaloSD.h:133
SensitiveDetectorCatalog
Definition: SensitiveDetectorCatalog.h:10
CaloSD::processHit
void processHit(const G4Step *step)
Definition: CaloSD.h:110
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
CaloSD::posGlobal
G4ThreeVector posGlobal
Definition: CaloSD.h:135
CaloSD::currentID
CaloHitID currentID
Definition: CaloSD.h:139
ALCARECOTkAlJpsiMuMu_cff.charge
charge
Definition: ALCARECOTkAlJpsiMuMu_cff.py:47
CaloG4Hit::getEM
double getEM() const
Definition: CaloG4Hit.h:55
edm::ParameterSet
Definition: ParameterSet.h:47
a
double a
Definition: hdecay.h:119
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
ZdcSD::ZdcSD
ZdcSD(const std::string &, const SensitiveDetectorCatalog &, edm::ParameterSet const &, const SimTrackManager *)
Definition: ZdcSD.cc:29
CaloSD::incidentEnergy
float incidentEnergy
Definition: CaloSD.h:136
GeV
const double GeV
Definition: MathUtil.h:16
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
CaloSD::resetForNewPrimary
void resetForNewPrimary(const G4Step *)
Definition: CaloSD.cc:670
IdealGeometryRecord.h
TrackInformation.h
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
generator_cfi.scheme
scheme
Definition: generator_cfi.py:22
alignCSCRings.r
r
Definition: alignCSCRings.py:93
DDAxes::phi
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
ZdcSD::setDetUnitId
uint32_t setDetUnitId(const G4Step *step) override
Definition: ZdcSD.cc:301
CaloSD::setParameterized
void setParameterized(bool val)
Definition: CaloSD.h:107
CaloG4Hit::getIncidentEnergy
double getIncidentEnergy() const
Definition: CaloG4Hit.h:61
CaloSD::entranceLocal
G4ThreeVector entranceLocal
Definition: CaloSD.h:134
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
ZdcSD::zdcHitEnergyCut
double zdcHitEnergyCut
Definition: ZdcSD.h:34
ztail.d
d
Definition: ztail.py:151
pi
const Double_t pi
Definition: trackSplitPlot.h:36
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
ZdcSD::initRun
void initRun() override
Definition: ZdcSD.cc:63
ParameterSet.h
PbPb_ZMuSkimMuonDPG_cff.particleType
particleType
Definition: PbPb_ZMuSkimMuonDPG_cff.py:27
ZdcSD::useShowerLibrary
bool useShowerLibrary
Definition: ZdcSD.h:32
ZdcSD::showerLibrary
std::unique_ptr< ZdcShowerLibrary > showerLibrary
Definition: ZdcSD.h:36
CaloG4Hit::getHadr
double getHadr() const
Definition: CaloG4Hit.h:58
ZdcSD::getEnergyDeposit
double getEnergyDeposit(const G4Step *) override
Definition: ZdcSD.cc:119
CaloSD
Definition: CaloSD.h:39
ZdcSD::getFromLibrary
bool getFromLibrary(const G4Step *) override
Definition: ZdcSD.cc:65
Generated for CMSSW Reference Manual by doxygen 1.8.16