Skip to content

Commit b22a038

Browse files
author
Sunanda
committed
Restore the Run1 code for ZDC hit producion usng a library paramerization which can be switched on through parameter setting
1 parent b8e3f2c commit b22a038

File tree

2 files changed

+253
-9
lines changed

2 files changed

+253
-9
lines changed

‎SimG4CMS/Forward/interface/ZdcSD.h‎

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,8 @@
99
#include "SimG4CMS/Calo/interface/CaloSD.h"
1010
#include "SimG4CMS/Forward/interface/ZdcShowerLibrary.h"
1111
#include "SimG4CMS/Forward/interface/ZdcNumberingScheme.h"
12+
#include "FWCore/Framework/interface/Frameworkfwd.h"
13+
#include "FWCore/ParameterSet/interface/ParameterSetfwd.h"
1214

1315
class ZdcSD : public CaloSD {
1416
public:
@@ -18,7 +20,10 @@ class ZdcSD : public CaloSD {
1820
bool ProcessHits(G4Step *step, G4TouchableHistory *tHistory) override;
1921
uint32_t setDetUnitId(const G4Step *step) override;
2022

21-
// bool getFromLibrary(const G4Step * step) override;
23+
protected:
24+
double getEnergyDeposit(const G4Step *) override;
25+
bool getFromLibrary(const G4Step *) override;
26+
void initRun() override;
2227

2328
double calculateCherenkovDeposit(const G4Step *);
2429
double calculateMeanNumberOfPhotons(int, double, double);
@@ -31,13 +36,11 @@ class ZdcSD : public CaloSD {
3136
double evaluateFunction(const std::vector<double> &, const std::vector<double> &, double);
3237
double linearInterpolation(double, double, double, double, double);
3338

34-
protected:
35-
void initRun() override;
39+
int setTrackID(const G4Step *step) override;
3640

3741
private:
3842
int verbosity;
3943
bool useShowerLibrary, useShowerHits;
40-
int setTrackID(const G4Step *step) override;
4144
double thFibDir;
4245
double zdcHitEnergyCut;
4346
std::unique_ptr<ZdcShowerLibrary> showerLibrary;

‎SimG4CMS/Forward/src/ZdcSD.cc‎

Lines changed: 246 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -7,11 +7,13 @@
77
#include <memory>
88

99
#include "SimG4CMS/Forward/interface/ZdcSD.h"
10+
#include "SimG4CMS/Forward/interface/ForwardName.h"
1011
#include "FWCore/MessageLogger/interface/MessageLogger.h"
1112
#include "FWCore/ParameterSet/interface/ParameterSet.h"
1213
#include "Geometry/Records/interface/IdealGeometryRecord.h"
1314
#include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
1415
#include "SimG4Core/Notification/interface/TrackInformation.h"
16+
#include "SimG4CMS/Forward/interface/ZdcNumberingScheme.h"
1517

1618
#include "G4SDManager.hh"
1719
#include "G4Step.hh"
@@ -27,7 +29,7 @@
2729
#include "G4Poisson.hh"
2830
#include "G4TwoVector.hh"
2931

30-
//#define EDM_ML_DEBUG
32+
#define EDM_ML_DEBUG
3133

3234
ZdcSD::ZdcSD(const std::string& name,
3335
const SensitiveDetectorCatalog& clg,
@@ -37,7 +39,8 @@ ZdcSD::ZdcSD(const std::string& name,
3739
edm::ParameterSet m_ZdcSD = p.getParameter<edm::ParameterSet>("ZdcSD");
3840
useShowerLibrary = m_ZdcSD.getParameter<bool>("UseShowerLibrary");
3941
useShowerHits = m_ZdcSD.getParameter<bool>("UseShowerHits");
40-
zdcHitEnergyCut = m_ZdcSD.getParameter<double>("ZdcHitEnergyCut") * GeV;
42+
zdcHitEnergyCut = m_ZdcSD.getParameter<double>("ZdcHitEnergyCut") * CLHEP::GeV;
43+
thFibDir = m_ZdcSD.getParameter<double>("FiberDirection");
4144
verbosity = m_ZdcSD.getParameter<int>("Verbosity");
4245
int verbn = verbosity / 10;
4346
verbosity %= 10;
@@ -73,10 +76,9 @@ void ZdcSD::initRun() {
7376
bool ZdcSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
7477
NaNTrap(aStep);
7578

76-
/*
77-
if (useShowerLibrary)
79+
if (useShowerLibrary)
7880
getFromLibrary(aStep);
79-
*/
81+
8082
#ifdef EDM_ML_DEBUG
8183
edm::LogVerbatim("ZdcSD") << "ZdcSD::" << GetName() << " ID= " << aStep->GetTrack()->GetTrackID()
8284
<< " prID= " << aStep->GetTrack()->GetParentID()
@@ -114,6 +116,245 @@ bool ZdcSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
114116
return true;
115117
}
116118

119+
bool ZdcSD::getFromLibrary(const G4Step* aStep) {
120+
bool ok = true;
121+
122+
auto const preStepPoint = aStep->GetPreStepPoint();
123+
124+
double etrack = preStepPoint->GetKineticEnergy();
125+
int primaryID = setTrackID(aStep);
126+
127+
hits.clear();
128+
129+
// Reset entry point for new primary
130+
resetForNewPrimary(aStep);
131+
132+
if (etrack >= zdcHitEnergyCut) {
133+
// create hits only if above threshold
134+
135+
#ifdef EDM_ML_DEBUG
136+
auto const theTrack = aStep->GetTrack();
137+
edm::LogVerbatim("ForwardSim") << "----------------New track------------------------------\n"
138+
<< "Incident EnergyTrack: " << etrack << " MeV \n"
139+
<< "Zdc Cut Energy for Hits: " << zdcHitEnergyCut << " MeV \n"
140+
<< "ZdcSD::getFromLibrary " << hits.size() << " hits for " << GetName() << " of "
141+
<< primaryID << " with " << theTrack->GetDefinition()->GetParticleName() << " of "
142+
<< etrack << " MeV\n";
143+
#endif
144+
hits.swap(showerLibrary.get()->getHits(aStep, ok));
145+
}
146+
147+
incidentEnergy = etrack;
148+
entrancePoint = preStepPoint->GetPosition();
149+
for (unsigned int i = 0; i < hits.size(); i++) {
150+
posGlobal = hits[i].position;
151+
entranceLocal = hits[i].entryLocal;
152+
double time = hits[i].time;
153+
unsigned int unitID = hits[i].detID;
154+
edepositHAD = hits[i].DeHad;
155+
edepositEM = hits[i].DeEM;
156+
currentID[0].setID(unitID, time, primaryID, 0);
157+
processHit(aStep);
158+
159+
#ifdef EDM_ML_DEBUG
160+
edm::LogVerbatim("ForwardSim") << "ZdcSD: Final Hit number:" << i << "-->"
161+
<< "New HitID: " << currentHit[0]->getUnitID()
162+
<< " New Hit trackID: " << currentHit[0]->getTrackID()
163+
<< " New EM Energy: " << currentHit[0]->getEM() / CLHEP::GeV
164+
<< " New HAD Energy: " << currentHit[0]->getHadr() / CLHEP::GeV
165+
<< " New HitEntryPoint: " << currentHit[0]->getEntryLocal()
166+
<< " New IncidentEnergy: " << currentHit[0]->getIncidentEnergy() / CLHEP::GeV
167+
<< " New HitPosition: " << posGlobal;
168+
#endif
169+
}
170+
return ok;
171+
}
172+
173+
double ZdcSD::getEnergyDeposit(const G4Step* aStep) {
174+
double NCherPhot = 0.;
175+
176+
// preStepPoint information
177+
G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
178+
G4VPhysicalVolume* currentPV = preStepPoint->GetPhysicalVolume();
179+
std::string nameVolume = ForwardName::getName(currentPV->GetName());
180+
181+
const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
182+
const G4ThreeVector& hit_mom = preStepPoint->GetMomentumDirection();
183+
G4double stepL = aStep->GetStepLength() / cm;
184+
G4double beta = preStepPoint->GetBeta();
185+
G4double charge = preStepPoint->GetCharge();
186+
187+
// theTrack information
188+
G4Track* theTrack = aStep->GetTrack();
189+
G4String particleType = theTrack->GetDefinition()->GetParticleName();
190+
G4ThreeVector localPoint = theTrack->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
191+
192+
#ifdef EDM_ML_DEBUG
193+
const G4ThreeVector& vert_mom = theTrack->GetVertexMomentumDirection();
194+
195+
// calculations
196+
float costheta =
197+
vert_mom.z() / sqrt(vert_mom.x() * vert_mom.x() + vert_mom.y() * vert_mom.y() + vert_mom.z() * vert_mom.z());
198+
float theta = std::acos(std::min(std::max(costheta, -1.f), 1.f));
199+
float eta = -std::log(std::tan(theta * 0.5f));
200+
float phi = -100.;
201+
if (vert_mom.x() != 0)
202+
phi = std::atan2(vert_mom.y(), vert_mom.x());
203+
if (phi < 0.)
204+
phi += twopi;
205+
206+
// Get the total energy deposit
207+
double stepE = aStep->GetTotalEnergyDeposit();
208+
209+
// postStepPoint information
210+
G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
211+
G4VPhysicalVolume* postPV = postStepPoint->GetPhysicalVolume();
212+
std::string postnameVolume = ForwardName::getName(postPV->GetName());
213+
edm::LogVerbatim("ForwardSim") << "ZdcSD:: getEnergyDeposit: \n"
214+
<< " preStepPoint: " << nameVolume << "," << stepL << "," << stepE << "," << beta
215+
<< "," << charge << "\n"
216+
<< " postStepPoint: " << postnameVolume << "," << costheta << "," << theta << ","
217+
<< eta << "," << phi << "," << particleType << " id= " << theTrack->GetTrackID()
218+
<< " Etot(GeV)= " << theTrack->GetTotalEnergy() / CLHEP::GeV;
219+
#endif
220+
const double bThreshold = 0.67;
221+
if ((beta > bThreshold) && (charge != 0) && (nameVolume == "ZDC_EMFiber" || nameVolume == "ZDC_HadFiber")) {
222+
#ifdef EDM_ML_DEBUG
223+
edm::LogVerbatim("ForwardSim") << "ZdcSD:: getEnergyDeposit: pass ";
224+
#endif
225+
const float nMedium = 1.4925;
226+
// float photEnSpectrDL = 10714.285714;
227+
// photEnSpectrDL = (1./400.nm-1./700.nm)*10000000.cm/nm; /* cm-1 */
228+
229+
const float photEnSpectrDE = 1.24;
230+
// E = 2pi*(1./137.)*(eV*cm/370.)/lambda = 12.389184*(eV*cm)/lambda
231+
// Emax = 12.389184*(eV*cm)/400nm*10-7cm/nm = 3.01 eV
232+
// Emin = 12.389184*(eV*cm)/700nm*10-7cm/nm = 1.77 eV
233+
// delE = Emax - Emin = 1.24 eV
234+
235+
const float effPMTandTransport = 0.15;
236+
237+
// Check these values
238+
const float thFullRefl = 23.;
239+
float thFullReflRad = thFullRefl * pi / 180.;
240+
241+
float thFibDirRad = thFibDir * pi / 180.;
242+
243+
// at which theta the point is located:
244+
// float th1 = hitPoint.theta();
245+
246+
// theta of charged particle in LabRF(hit momentum direction):
247+
float costh = hit_mom.z() / sqrt(hit_mom.x() * hit_mom.x() + hit_mom.y() * hit_mom.y() + hit_mom.z() * hit_mom.z());
248+
float th = acos(std::min(std::max(costh, -1.f), 1.f));
249+
// just in case (can do both standard ranges of phi):
250+
if (th < 0.)
251+
th += CLHEP::twopi;
252+
253+
// theta of cone with Cherenkov photons w.r.t.direction of charged part.:
254+
float costhcher = 1. / (nMedium * beta);
255+
float thcher = acos(std::min(std::max(costhcher, -1.f), 1.f));
256+
257+
// diff thetas of charged part. and quartz direction in LabRF:
258+
float DelFibPart = std::abs(th - thFibDirRad);
259+
260+
// define real distances:
261+
float d = std::abs(std::tan(th) - std::tan(thFibDirRad));
262+
263+
float a = std::tan(thFibDirRad) + std::tan(std::abs(thFibDirRad - thFullReflRad));
264+
float r = std::tan(th) + std::tan(std::abs(th - thcher));
265+
266+
// define losses d_qz in cone of full reflection inside quartz direction
267+
float d_qz = -1;
268+
#ifdef EDM_ML_DEBUG
269+
float variant = -1;
270+
#endif
271+
// if (d > (r+a))
272+
if (DelFibPart > (thFullReflRad + thcher)) {
273+
#ifdef EDM_ML_DEBUG
274+
variant = 0.;
275+
#endif
276+
d_qz = 0.;
277+
} else {
278+
// if ((DelFibPart + thcher) < thFullReflRad ) [(d+r) < a]
279+
if ((th + thcher) < (thFibDirRad + thFullReflRad) && (th - thcher) > (thFibDirRad - thFullReflRad)) {
280+
#ifdef EDM_ML_DEBUG
281+
variant = 1.;
282+
#endif
283+
d_qz = 1.;
284+
} else {
285+
// if ((thcher - DelFibPart ) > thFullReflRad ) [(r-d) > a]
286+
if ((thFibDirRad + thFullReflRad) < (th + thcher) && (thFibDirRad - thFullReflRad) > (th - thcher)) {
287+
#ifdef EDM_ML_DEBUG
288+
variant = 2.;
289+
#endif
290+
d_qz = 0.;
291+
} else {
292+
#ifdef EDM_ML_DEBUG
293+
variant = 3.; // d_qz is calculated below
294+
#endif
295+
// use crossed length of circles(cone projection) - dC1/dC2 :
296+
float arg_arcos = 0.;
297+
float tan_arcos = 2. * a * d;
298+
if (tan_arcos != 0.)
299+
arg_arcos = (r * r - a * a - d * d) / tan_arcos;
300+
arg_arcos = std::abs(arg_arcos);
301+
float th_arcos = acos(std::min(std::max(arg_arcos, -1.f), 1.f));
302+
d_qz = th_arcos / CLHEP::twopi;
303+
d_qz = std::abs(d_qz);
304+
#ifdef EDM_ML_DEBUG
305+
edm::LogVerbatim("ForwardSim") << " d_qz: " << r << "," << a << "," << d << " " << tan_arcos << " "
306+
<< arg_arcos;
307+
edm::LogVerbatim("ForwardSim") << "," << arg_arcos;
308+
edm::LogVerbatim("ForwardSim") << " " << d_qz;
309+
edm::LogVerbatim("ForwardSim") << " " << th_arcos;
310+
edm::LogVerbatim("ForwardSim") << "," << d_qz;
311+
#endif
312+
}
313+
}
314+
}
315+
double meanNCherPhot = 0.;
316+
int poissNCherPhot = 0;
317+
if (d_qz > 0) {
318+
meanNCherPhot = 370. * charge * charge * (1. - 1. / (nMedium * nMedium * beta * beta)) * photEnSpectrDE * stepL;
319+
320+
poissNCherPhot = std::max((int)G4Poisson(meanNCherPhot), 0);
321+
NCherPhot = poissNCherPhot * effPMTandTransport * d_qz;
322+
}
323+
324+
#ifdef EDM_ML_DEBUG
325+
edm::LogVerbatim("ForwardSim") << "ZdcSD:: getEnergyDeposit: gED: " << stepE << "," << costh << "," << th << ","
326+
<< costhcher << "," << thcher << "," << DelFibPart << "," << d << "," << a << ","
327+
<< r << "," << hitPoint << "," << hit_mom << "," << vert_mom << "," << localPoint
328+
<< "," << charge << "," << beta << "," << stepL << "," << d_qz << "," << variant
329+
<< "," << meanNCherPhot << "," << poissNCherPhot << "," << NCherPhot;
330+
#endif
331+
// --constants-----------------
332+
// << "," << photEnSpectrDE
333+
// << "," << nMedium
334+
// << "," << bThreshold
335+
// << "," << thFibDirRad
336+
// << "," << thFullReflRad
337+
// << "," << effPMTandTransport
338+
// --other variables-----------
339+
// << "," << curprocess
340+
// << "," << nameProcess
341+
// << "," << name
342+
// << "," << rad
343+
// << "," << mat
344+
345+
} else {
346+
// determine failure mode: beta, charge, and/or nameVolume
347+
if (beta <= bThreshold)
348+
edm::LogVerbatim("ForwardSim") << "ZdcSD:: getEnergyDeposit: fail beta=" << beta;
349+
if (charge == 0)
350+
edm::LogVerbatim("ForwardSim") << "ZdcSD:: getEnergyDeposit: fail charge=0";
351+
if (!(nameVolume == "ZDC_EMFiber" || nameVolume == "ZDC_HadFiber"))
352+
edm::LogVerbatim("ForwardSim") << "ZdcSD:: getEnergyDeposit: fail nv=" << nameVolume;
353+
}
354+
355+
return NCherPhot;
356+
}
357+
117358
///////////////////////////////////////
118359
// Functions added by Oliver Suranyi //
119360
///////////////////////////////////////

0 commit comments

Comments
 (0)