CMS 3D CMS Logo

List of all members | Public Member Functions | Protected Member Functions | Private Member Functions | Private Attributes
CastorSD Class Reference

#include <SimG4CMS/Forward/interface/CastorSD.h>

Inheritance diagram for CastorSD:
CaloSD SensitiveCaloDetector Observer< const BeginOfRun *> Observer< const BeginOfEvent *> Observer< const BeginOfTrack *> Observer< const EndOfTrack *> Observer< const EndOfEvent *> SensitiveDetector

Public Member Functions

 CastorSD (const std::string &, const SensitiveDetectorCatalog &clg, edm::ParameterSet const &, const SimTrackManager *)
 
uint32_t setDetUnitId (const G4Step *step) override
 
void setNumberingScheme (CastorNumberingScheme *scheme)
 
 ~CastorSD () override
 
- Public Member Functions inherited from CaloSD
 CaloSD (const std::string &aSDname, const SensitiveDetectorCatalog &clg, edm::ParameterSet const &p, const SimTrackManager *, float timeSlice=1., bool ignoreTkID=false)
 
void clear () override
 
void clearHits () override
 
void DrawAll () override
 
void EndOfEvent (G4HCofThisEvent *eventHC) override
 
void fillHits (edm::PCaloHitContainer &, const std::string &) override
 
void Initialize (G4HCofThisEvent *HCE) override
 
bool isItFineCalo (const G4VTouchable *touch)
 
void PrintAll () override
 
G4bool ProcessHits (G4Step *step, G4TouchableHistory *) override
 
bool ProcessHits (G4GFlashSpot *aSpot, G4TouchableHistory *) override
 
void reset () override
 
 ~CaloSD () override
 
- Public Member Functions inherited from SensitiveCaloDetector
 SensitiveCaloDetector (const std::string &iname, const SensitiveDetectorCatalog &clg)
 
- Public Member Functions inherited from SensitiveDetector
void EndOfEvent (G4HCofThisEvent *eventHC) override
 
const std::vector< std::string > & getNames () const
 
void Initialize (G4HCofThisEvent *eventHC) override
 
bool isCaloSD () const
 
 SensitiveDetector (const std::string &iname, const SensitiveDetectorCatalog &, bool calo)
 
 ~SensitiveDetector () override
 
- Public Member Functions inherited from Observer< const BeginOfRun *>
 Observer ()
 
void slotForUpdate (const BeginOfRun * iT)
 
virtual ~Observer ()
 
- Public Member Functions inherited from Observer< const BeginOfEvent *>
 Observer ()
 
void slotForUpdate (const BeginOfEvent * iT)
 
virtual ~Observer ()
 
- Public Member Functions inherited from Observer< const BeginOfTrack *>
 Observer ()
 
void slotForUpdate (const BeginOfTrack * iT)
 
virtual ~Observer ()
 
- Public Member Functions inherited from Observer< const EndOfTrack *>
 Observer ()
 
void slotForUpdate (const EndOfTrack * iT)
 
virtual ~Observer ()
 
- Public Member Functions inherited from Observer< const EndOfEvent *>
 Observer ()
 
void slotForUpdate (const EndOfEvent * iT)
 
virtual ~Observer ()
 

Protected Member Functions

double getEnergyDeposit (const G4Step *) override
 
bool getFromLibrary (const G4Step *) override
 
- Protected Member Functions inherited from CaloSD
bool checkHit ()
 
CaloG4HitcreateNewHit (const G4Step *, const G4Track *)
 
virtual void endEvent ()
 
virtual double EnergyCorrected (const G4Step &step, const G4Track *)
 
virtual bool filterHit (CaloG4Hit *, double)
 
unsigned int findBoundaryCrossingParent (const G4Track *track, bool markParentAsSaveable=true)
 
double getAttenuation (const G4Step *aStep, double birk1, double birk2, double birk3) const
 
virtual uint16_t getDepth (const G4Step *)
 
int getNumberOfHits ()
 
double getResponseWt (const G4Track *)
 
virtual int getTrackID (const G4Track *)
 
bool hitExists (const G4Step *)
 
void ignoreRejection ()
 
virtual void initEvent (const BeginOfEvent *)
 
virtual void initRun ()
 
void printDetectorLevels (const G4VTouchable *) const
 
void processHit (const G4Step *step)
 
void resetForNewPrimary (const G4Step *)
 
void setNumberCheckedHits (int val)
 
void setParameterized (bool val)
 
G4ThreeVector setToGlobal (const G4ThreeVector &, const G4VTouchable *) const
 
G4ThreeVector setToLocal (const G4ThreeVector &, const G4VTouchable *) const
 
virtual int setTrackID (const G4Step *)
 
void setUseMap (bool val)
 
std::string shortreprID (const CaloHitID &ID)
 
std::string shortreprID (const CaloG4Hit *hit)
 
void update (const BeginOfRun *) override
 This routine will be called when the appropriate signal arrives. More...
 
void update (const BeginOfEvent *) override
 This routine will be called when the appropriate signal arrives. More...
 
void update (const BeginOfTrack *trk) override
 This routine will be called when the appropriate signal arrives. More...
 
void update (const EndOfTrack *trk) override
 This routine will be called when the appropriate signal arrives. More...
 
void update (const ::EndOfEvent *) override
 
void updateHit (CaloG4Hit *)
 
- Protected Member Functions inherited from SensitiveDetector
TrackInformationcmsTrackInformation (const G4Track *aTrack)
 
Local3DPoint ConvertToLocal3DPoint (const G4ThreeVector &point) const
 
Local3DPoint FinalStepPosition (const G4Step *step, coordinates) const
 
Local3DPoint InitialStepPosition (const G4Step *step, coordinates) const
 
Local3DPoint LocalPostStepPosition (const G4Step *step) const
 
Local3DPoint LocalPreStepPosition (const G4Step *step) const
 
void NaNTrap (const G4Step *step) const
 
void setNames (const std::vector< std::string > &)
 
- Protected Member Functions inherited from Observer< const EndOfEvent *>
virtual void update (const EndOfEvent *)=0
 This routine will be called when the appropriate signal arrives. More...
 

Private Member Functions

uint32_t rotateUnitID (uint32_t, const G4Track *, const CastorShowerEvent &)
 

Private Attributes

double energyThresholdSL
 
G4LogicalVolume * lvC3EF
 
G4LogicalVolume * lvC3HF
 
G4LogicalVolume * lvC4EF
 
G4LogicalVolume * lvC4HF
 
G4LogicalVolume * lvCAST
 
double non_compensation_factor
 
CastorNumberingSchemenumberingScheme
 
CastorShowerLibraryshowerLibrary
 
bool useShowerLibrary
 

Additional Inherited Members

- Protected Types inherited from SensitiveDetector
enum  coordinates { WorldCoordinates, LocalCoordinates }
 
- Static Protected Member Functions inherited from CaloSD
static std::string printableDecayChain (const std::vector< unsigned int > &decayChain)
 
- Protected Attributes inherited from CaloSD
CaloG4HitcurrentHit
 
CaloHitID currentID
 
float edepositEM
 
float edepositHAD
 
double eminHit
 
double energyCut
 
G4ThreeVector entranceLocal
 
G4ThreeVector entrancePoint
 
bool forceSave
 
float incidentEnergy
 
double kmaxIon
 
double kmaxNeutron
 
double kmaxProton
 
G4ThreeVector posGlobal
 
CaloHitID previousID
 
bool suppressHeavy
 
double tmaxHit
 

Detailed Description

Description: Stores hits of Castor in appropriate container

Usage: Used in sensitive detector builder

Definition at line 33 of file CastorSD.h.

Constructor & Destructor Documentation

◆ CastorSD()

CastorSD::CastorSD ( const std::string &  name,
const SensitiveDetectorCatalog clg,
edm::ParameterSet const &  p,
const SimTrackManager manager 
)

Definition at line 31 of file CastorSD.cc.

References g4SimHits_cfi::CastorShowerLibrary, energyThresholdSL, edm::ParameterSet::getParameter(), lvC3EF, lvC3HF, lvC4EF, lvC4HF, lvCAST, Skims_PA_cff::name, non_compensation_factor, AlCaHLTBitMon_ParallelJobs::p, setNumberingScheme(), CaloSD::setParameterized(), showerLibrary, and useShowerLibrary.

35  : CaloSD(name, clg, p, manager),
36  numberingScheme(nullptr),
37  lvC3EF(nullptr),
38  lvC3HF(nullptr),
39  lvC4EF(nullptr),
40  lvC4HF(nullptr),
41  lvCAST(nullptr) {
42  edm::ParameterSet m_CastorSD = p.getParameter<edm::ParameterSet>("CastorSD");
43  useShowerLibrary = m_CastorSD.getParameter<bool>("useShowerLibrary");
44  energyThresholdSL = m_CastorSD.getParameter<double>("minEnergyInGeVforUsingSLibrary");
45  energyThresholdSL = energyThresholdSL * CLHEP::GeV; // Convert GeV => MeV
46 
47  non_compensation_factor = m_CastorSD.getParameter<double>("nonCompensationFactor");
48 
49  if (useShowerLibrary) {
51  setParameterized(true);
52  }
54 
55  edm::LogVerbatim("ForwardSim") << "********************************************************\n"
56  << "* Constructing a CastorSD with name " << GetName() << "\n"
57  << "* Using Castor Shower Library: " << useShowerLibrary << "\n"
58  << "********************************************************";
59 
60  const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
61  for (auto lv : *lvs) {
62  if (strcmp((lv->GetName()).c_str(), "C3EF") == 0) {
63  lvC3EF = lv;
64  }
65  if (strcmp((lv->GetName()).c_str(), "C3HF") == 0) {
66  lvC3HF = lv;
67  }
68  if (strcmp((lv->GetName()).c_str(), "C4EF") == 0) {
69  lvC4EF = lv;
70  }
71  if (strcmp((lv->GetName()).c_str(), "C4HF") == 0) {
72  lvC4HF = lv;
73  }
74  if (strcmp((lv->GetName()).c_str(), "CAST") == 0) {
75  lvCAST = lv;
76  }
77  if (lvC3EF != nullptr && lvC3HF != nullptr && lvC4EF != nullptr && lvC4HF != nullptr && lvCAST != nullptr) {
78  break;
79  }
80  }
81  edm::LogVerbatim("ForwardSim") << "CastorSD:: LogicalVolume pointers\n"
82  << lvC3EF << " for C3EF; " << lvC3HF << " for C3HF; " << lvC4EF << " for C4EF; "
83  << lvC4HF << " for C4HF; " << lvCAST << " for CAST. ";
84 }
Log< level::Info, true > LogVerbatim
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
double non_compensation_factor
Definition: CastorSD.h:53
CastorNumberingScheme * numberingScheme
Definition: CastorSD.h:46
bool useShowerLibrary
Definition: CastorSD.h:51
G4LogicalVolume * lvC4HF
Definition: CastorSD.h:48
G4LogicalVolume * lvC3HF
Definition: CastorSD.h:48
double energyThresholdSL
Definition: CastorSD.h:52
G4LogicalVolume * lvC4EF
Definition: CastorSD.h:48
G4LogicalVolume * lvCAST
Definition: CastorSD.h:49
CastorShowerLibrary * showerLibrary
Definition: CastorSD.h:47
void setNumberingScheme(CastorNumberingScheme *scheme)
Definition: CastorSD.cc:289
G4LogicalVolume * lvC3EF
Definition: CastorSD.h:48
CaloSD(const std::string &aSDname, const SensitiveDetectorCatalog &clg, edm::ParameterSet const &p, const SimTrackManager *, float timeSlice=1., bool ignoreTkID=false)
Definition: CaloSD.cc:33
void setParameterized(bool val)
Definition: CaloSD.h:110

◆ ~CastorSD()

CastorSD::~CastorSD ( )
override

Definition at line 88 of file CastorSD.cc.

References showerLibrary.

88 { delete showerLibrary; }
CastorShowerLibrary * showerLibrary
Definition: CastorSD.h:47

Member Function Documentation

◆ getEnergyDeposit()

double CastorSD::getEnergyDeposit ( const G4Step *  aStep)
overrideprotectedvirtual

Reimplemented from CaloSD.

Definition at line 92 of file CastorSD.cc.

References a, funct::abs(), HLT_2022v12_cff::beta, ALCARECOTkAlJpsiMuMu_cff::charge, ztail::d, TrackInformation::getCastorHitPID(), TrackInformation::hasCastorHit(), G4TrackToParticleID::isStableHadronIon(), lvC3EF, lvC3HF, lvC4EF, lvC4HF, SiStripPI::max, SiStripPI::min, non_compensation_factor, pi, alignCSCRings::r, L1EGammaClusterEmuProducer_cfi::scale, TrackInformation::setCastorHitPID(), mathSSE::sqrt(), and funct::tan().

92  {
93  double NCherPhot = 0.;
94 
95  // Get theTrack
96  auto const theTrack = aStep->GetTrack();
97 
98  // preStepPoint information *********************************************
99  auto const preStepPoint = aStep->GetPreStepPoint();
100  auto const currentPV = preStepPoint->GetPhysicalVolume();
101  auto const currentLV = currentPV->GetLogicalVolume();
102 
103 #ifdef EDM_ML_DEBUG
104  edm::LogVerbatim("ForwardSim") << "CastorSD::getEnergyDeposit:"
105  << "\n CurrentStepNumber , TrackID , ParentID, Particle , VertexPosition ,"
106  << " LogicalVolumeAtVertex , PV, Time"
107  << "\n TRACKINFO: " << theTrack->GetCurrentStepNumber() << " , "
108  << theTrack->GetTrackID() << " , " << theTrack->GetParentID() << " , "
109  << theTrack->GetDefinition()->GetParticleName() << " , "
110  << theTrack->GetVertexPosition() << " , "
111  << theTrack->GetLogicalVolumeAtVertex()->GetName() << " , " << currentPV->GetName()
112  << " , " << theTrack->GetGlobalTime();
113 #endif
114 
115  // if particle moves from interaction point or "backwards (halo)
116  const G4ThreeVector& hit_mom = preStepPoint->GetMomentumDirection();
117  const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
118  double zint = hitPoint.z();
119 
120  // Check if theTrack is a muon (if so, DO NOT use Shower Library)
121  G4int parCode = theTrack->GetDefinition()->GetPDGEncoding();
122  bool isHad = G4TrackToParticleID::isStableHadronIon(theTrack);
123 
124  double meanNCherPhot = 0;
125  G4double charge = preStepPoint->GetCharge();
126  // VI: no Cerenkov light from neutrals
127  if (0.0 == charge) {
128  return meanNCherPhot;
129  }
130 
131  G4double beta = preStepPoint->GetBeta();
132  const double bThreshold = 0.67;
133  // VI: no Cerenkov light from non-relativistic particles
134  if (beta < bThreshold) {
135  return meanNCherPhot;
136  }
137 
138  // remember primary particle hitting the CASTOR detector
139  TrackInformationExtractor TIextractor;
140  TrackInformation& trkInfo = TIextractor(theTrack);
141  if (!trkInfo.hasCastorHit()) {
142  trkInfo.setCastorHitPID(parCode);
143  }
144  G4double stepl = aStep->GetStepLength() / CLHEP::cm;
145 
146  //int castorHitPID = trkInfo.hasCastorHit() ? std::abs(trkInfo.getCastorHitPID())
147  // : std::abs(parCode);
148 
149  // *************************************************
150  // take into account light collection curve for plate
151  // double weight = curve_Castor(nameVolume, preStepPoint);
152  // double edep = aStep->GetTotalEnergyDeposit() * weight;
153  // *************************************************
154 
155  // *************************************************
156  /* comments for sensitive volumes:
157  C001 ...-... CP06,CBU1 ...-...CALM --- > fibres and bundle
158  for first release of CASTOR
159  CASF --- > quartz plate for first and second releases of CASTOR
160  GF2Q, GFNQ, GR2Q, GRNQ
161  for tests with my own test geometry of HF (on ask of Gavrilov)
162  C3TF, C4TF - for third release of CASTOR
163  */
164 #ifdef EDM_ML_DEBUG
165  edm::LogVerbatim("ForwardSim") << "CastorSD::getEnergyDeposit: for ID=" << theTrack->GetTrackID()
166  << " LV: " << currentLV->GetName() << " isHad:" << isHad << " pdg=" << parCode
167  << " castorPID=" << trkInfo.getCastorHitPID() << " sl= " << stepl
168  << " Edep= " << aStep->GetTotalEnergyDeposit();
169 #endif
170  if (currentLV == lvC3EF || currentLV == lvC4EF || currentLV == lvC3HF || currentLV == lvC4HF) {
171  const double nMedium = 1.4925;
172  // double photEnSpectrDL = (1./400.nm-1./700.nm)*10000000.cm/nm; /* cm-1 */
173  // double photEnSpectrDL = 10714.285714;
174 
175  const double photEnSpectrDE = 1.24; /* see below */
176  /* E = 2pi*(1./137.)*(eV*cm/370.)/lambda = */
177  /* = 12.389184*(eV*cm)/lambda */
178  /* Emax = 12.389184*(eV*cm)/400nm*10-7cm/nm = 3.01 eV */
179  /* Emin = 12.389184*(eV*cm)/700nm*10-7cm/nm = 1.77 eV */
180  /* delE = Emax - Emin = 1.24 eV --> */
181  /* */
182  /* default for Castor nameVolume == "CASF" or (C3TF & C4TF) */
183 
184  const double thFullRefl = 23.; /* 23.dergee */
185  const double thFullReflRad = thFullRefl * pi / 180.;
186 
187  /* default for Castor nameVolume == "CASF" or (C3TF & C4TF) */
188  const double thFibDir = 45.; /* .dergee */
189  /* for test HF geometry volumes:
190  if(nameVolume == "GF2Q" || nameVolume == "GFNQ" ||
191  nameVolume == "GR2Q" || nameVolume == "GRNQ")
192  thFibDir = 0.0; // .dergee
193  */
194  const double thFibDirRad = thFibDir * pi / 180.;
195 
196  // theta of charged particle in LabRF(hit momentum direction):
197  double costh =
198  hit_mom.z() / sqrt(hit_mom.x() * hit_mom.x() + hit_mom.y() * hit_mom.y() + hit_mom.z() * hit_mom.z());
199  if (zint < 0)
200  costh = -costh;
201  double th = acos(std::min(std::max(costh, double(-1.)), double(1.)));
202 
203  // just in case (can do bot use):
204  if (th < 0.)
205  th += twopi;
206 
207  // theta of cone with Cherenkov photons w.r.t.direction of charged part.:
208  double costhcher = 1. / (nMedium * beta);
209  double thcher = acos(std::min(std::max(costhcher, double(-1.)), double(1.)));
210 
211  // diff thetas of charged part. and quartz direction in LabRF:
212  double DelFibPart = std::abs(th - thFibDirRad);
213 
214  // define real distances:
215  double d = std::abs(tan(th) - tan(thFibDirRad));
216 
217  double a = tan(thFibDirRad) + tan(std::abs(thFibDirRad - thFullReflRad));
218  double r = tan(th) + tan(std::abs(th - thcher));
219 
220  // define losses d_qz in cone of full reflection inside quartz direction
221  double d_qz;
222 #ifdef EDM_ML_DEBUG
223  int variant(0);
224 #endif
225  if (DelFibPart > (thFullReflRad + thcher)) {
226  d_qz = 0.;
227  } else {
228  if ((th + thcher) < (thFibDirRad + thFullReflRad) && (th - thcher) > (thFibDirRad - thFullReflRad)) {
229  d_qz = 1.;
230 #ifdef EDM_ML_DEBUG
231  variant = 1;
232 #endif
233  } else {
234  if ((thFibDirRad + thFullReflRad) < (th + thcher) && (thFibDirRad - thFullReflRad) > (th - thcher)) {
235  d_qz = 0.;
236 #ifdef EDM_ML_DEBUG
237  variant = 2;
238 #endif
239  } else {
240 #ifdef EDM_ML_DEBUG
241  variant = 3;
242 #endif
243  // use crossed length of circles(cone projection)
244  // dC1/dC2 :
245  double arg_arcos = 0.;
246  double tan_arcos = 2. * a * d;
247  if (tan_arcos != 0.)
248  arg_arcos = (r * r - a * a - d * d) / tan_arcos;
249  arg_arcos = std::abs(arg_arcos);
250  double th_arcos = acos(std::min(std::max(arg_arcos, -1.), 1.));
251  d_qz = std::abs(th_arcos / CLHEP::twopi);
252  }
253  }
254  }
255 #ifdef EDM_ML_DEBUG
256  edm::LogVerbatim("ForwardSim") << "variant" << variant;
257 #endif
258 
259  // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
260 
261  meanNCherPhot = 370. * charge * charge * (1. - 1. / (nMedium * nMedium * beta * beta)) * photEnSpectrDE * stepl;
262 
263  double scale = isHad ? non_compensation_factor : 1.0;
264 
265  int poissNCherPhot = std::max(0, (int)G4Poisson(meanNCherPhot * scale));
266 
267  const double effPMTandTransport = 0.19;
268  const double ReflPower = 0.1;
269  double proba = d_qz + (1 - d_qz) * ReflPower;
270  NCherPhot = poissNCherPhot * effPMTandTransport * proba * 0.307;
271 #ifdef EDM_ML_DEBUG
272  edm::LogVerbatim("ForwardSim") << " Nph= " << NCherPhot << " Np= " << poissNCherPhot
273  << " eff= " << effPMTandTransport << " pb= " << proba << " Nmean= " << meanNCherPhot
274  << " q=" << charge << " beta=" << beta << " nMedium= " << nMedium << " sl= " << stepl
275  << " Nde=" << photEnSpectrDE;
276 #endif
277  }
278  return NCherPhot;
279 }
Log< level::Info, true > LogVerbatim
int getCastorHitPID() const
bool hasCastorHit() const
double non_compensation_factor
Definition: CastorSD.h:53
G4LogicalVolume * lvC4HF
Definition: CastorSD.h:48
const Double_t pi
void setCastorHitPID(const int pid)
T sqrt(T t)
Definition: SSEVec.h:19
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
G4LogicalVolume * lvC3HF
Definition: CastorSD.h:48
d
Definition: ztail.py:151
static bool isStableHadronIon(const G4Track *)
G4LogicalVolume * lvC4EF
Definition: CastorSD.h:48
G4LogicalVolume * lvC3EF
Definition: CastorSD.h:48
double a
Definition: hdecay.h:119

◆ getFromLibrary()

bool CastorSD::getFromLibrary ( const G4Step *  aStep)
overrideprotectedvirtual

Reimplemented from CaloSD.

Definition at line 361 of file CastorSD.cc.

References funct::abs(), CaloSD::currentID, dot(), CaloSD::edepositEM, CaloSD::edepositHAD, energyThresholdSL, funct::false, CastorShowerLibrary::getShowerHits(), TrackInformation::hasCastorHit(), hfClusterShapes_cfi::hits, mps_fire::i, G4TrackToParticleID::isGammaElectronPositron(), G4TrackToParticleID::isStableHadronIon(), lvCAST, M_PI, non_compensation_factor, CaloSD::processHit(), dttmaxenums::R, CaloSD::resetForNewPrimary(), rotateUnitID(), L1EGammaClusterEmuProducer_cfi::scale, TrackInformation::setCastorHitPID(), CaloHitID::setID(), CaloSD::setTrackID(), showerLibrary, mathSSE::sqrt(), theta(), protons_cff::time, funct::true, and useShowerLibrary.

361  {
363  //
364  // Method to get hits from the Shower Library
365  //
366  // "updateHit" save the Hits to a CaloG4Hit container
367  //
369 
370  auto const theTrack = aStep->GetTrack();
371  auto parCode = theTrack->GetDefinition()->GetPDGEncoding();
372 
373  // remember primary particle hitting the CASTOR detector
374  TrackInformationExtractor TIextractor;
375  TrackInformation& trkInfo = TIextractor(theTrack);
376  if (!trkInfo.hasCastorHit()) {
377  trkInfo.setCastorHitPID(parCode);
378  }
379 
380  if (!useShowerLibrary) {
381  return false;
382  }
383 
384  // preStepPoint information *********************************************
385 
386  auto const preStepPoint = aStep->GetPreStepPoint();
387  auto const currentPV = preStepPoint->GetPhysicalVolume();
388  auto const currentLV = currentPV->GetLogicalVolume();
389 
390 #ifdef EDM_ML_DEBUG
391  edm::LogVerbatim("ForwardSim") << "CastorSD::getFromLibrary: for ID=" << theTrack->GetTrackID()
392  << " parentID= " << theTrack->GetParentID() << " "
393  << theTrack->GetDefinition()->GetParticleName() << " LV: " << currentLV->GetName()
394  << " PV: " << currentPV->GetName() << "\n eta= " << theTrack->GetPosition().eta()
395  << " phi= " << theTrack->GetPosition().phi()
396  << " z(cm)= " << theTrack->GetPosition().z() / CLHEP::cm
397  << " time(ns)= " << theTrack->GetGlobalTime()
398  << " E(GeV)= " << theTrack->GetTotalEnergy() / CLHEP::GeV;
399 
400 #endif
401 
402  const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
403  const G4ThreeVector& hit_mom = preStepPoint->GetMomentumDirection();
404  double zint = hitPoint.z();
405  double pz = hit_mom.z();
406 
407  // Check if theTrack moves backward
408  bool backward = (pz * zint < 0.) ? true : false;
409 
410  // Check that theTrack is above the energy threshold to use Shower Library
411  bool aboveThreshold = (theTrack->GetKineticEnergy() > energyThresholdSL) ? true : false;
412 
413  // Check if theTrack is a muon (if so, DO NOT use Shower Library)
415  bool isHad = G4TrackToParticleID::isStableHadronIon(theTrack);
416 
417  // angle condition
418  double theta_max = M_PI - 3.1305; // angle in radians corresponding to -5.2 eta
419  double R_mom = sqrt(hit_mom.x() * hit_mom.x() + hit_mom.y() * hit_mom.y());
420  double theta = atan2(R_mom, std::abs(pz));
421  bool angleok = (theta < theta_max) ? true : false;
422 
423  // OkToUse
424  double R = sqrt(hitPoint.x() * hitPoint.x() + hitPoint.y() * hitPoint.y());
425  bool dot = (zint < -14450. && R < 45.) ? true : false;
426  bool inRange = (zint < -14700. || R > 193.) ? false : true;
427  //bool OkToUse = (inRange && !dot) ? true : false;
428 
429  bool particleWithinShowerLibrary =
430  aboveThreshold && (isEM || isHad) && (!backward) && inRange && !dot && angleok && currentLV == lvCAST;
431 
432 #ifdef EDM_ML_DEBUG
433  edm::LogVerbatim("ForwardSim") << "CastorSD::getFromLibrary: ID= " << theTrack->GetTrackID() << " E>E0 "
434  << aboveThreshold << " isEM " << isEM << " isHad " << isHad << " backword " << backward
435  << " Ok " << (inRange && !dot) << " angle " << angleok
436  << " LV: " << currentLV->GetName() << " " << (currentLV == lvCAST) << " "
437  << particleWithinShowerLibrary << " Edep= " << aStep->GetTotalEnergyDeposit();
438 #endif
439 
440  // Use Castor shower library if energy is above threshold, is not a muon
441  // and is not moving backward
442  if (!particleWithinShowerLibrary) {
443  return false;
444  }
445 
446  // **** Call method to retrieve hits from the ShowerLibrary ****
447  // always kill primary
448  bool isKilled(true);
450 
451  int primaryID = setTrackID(aStep);
452 
453  // Reset entry point for new primary
454  resetForNewPrimary(aStep);
455 
456 #ifdef EDM_ML_DEBUG
457  edm::LogVerbatim("ForwardSim") << "CastorSD::getFromLibrary: " << hits.getNhit() << " hits for " << GetName()
458  << " from " << theTrack->GetDefinition()->GetParticleName() << " of "
459  << preStepPoint->GetKineticEnergy() / CLHEP::GeV << " GeV and trackID "
460  << theTrack->GetTrackID() << " isHad: " << isHad;
461 #endif
462 
463  // Scale to correct energy
464  double E_track = preStepPoint->GetTotalEnergy();
465  double E_SLhit = hits.getPrimE() * CLHEP::GeV;
466  double scale = E_track / E_SLhit;
467 
468  //Non compensation
469  if (isHad) {
470  scale *= non_compensation_factor; // if hadronic extend the scale with the non-compensation factor
471  }
472  // Loop over hits retrieved from the library
473  for (unsigned int i = 0; i < hits.getNhit(); ++i) {
474  // Get nPhotoElectrons and set edepositEM / edepositHAD accordingly
475  double nPhotoElectrons = hits.getNphotons(i) * scale;
476 
477  if (isEM) {
478  edepositEM = nPhotoElectrons;
479  edepositHAD = 0.;
480  } else {
481  edepositEM = 0.;
482  edepositHAD = nPhotoElectrons;
483  }
484 
485  // Get hit position and time
486  double time = hits.getTime(i);
487 
488  // Get hit detID
489  unsigned int unitID = hits.getDetID(i);
490 
491  // Make the detID "rotation" from one sector to another taking into account the
492  // sectors of the impinging particle (theTrack) and of the particle that produced
493  // the 'hits' retrieved from shower library
494  unsigned int rotatedUnitID = rotateUnitID(unitID, theTrack, hits);
495  currentID.setID(rotatedUnitID, time, primaryID, 0);
496  processHit(aStep);
497  } // End of loop over hits
498  return isKilled;
499 }
float edepositEM
Definition: CaloSD.h:140
Log< level::Info, true > LogVerbatim
bool hasCastorHit() const
double non_compensation_factor
Definition: CastorSD.h:53
bool useShowerLibrary
Definition: CastorSD.h:51
void processHit(const G4Step *step)
Definition: CaloSD.h:113
T dot(const Basic3DVector &v) const
Scalar product, or "dot" product, with a vector of same type.
void setCastorHitPID(const int pid)
CastorShowerEvent getShowerHits(const G4Step *, bool &)
float edepositHAD
Definition: CaloSD.h:140
T sqrt(T t)
Definition: SSEVec.h:19
void resetForNewPrimary(const G4Step *)
Definition: CaloSD.cc:644
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double energyThresholdSL
Definition: CastorSD.h:52
#define M_PI
static bool isStableHadronIon(const G4Track *)
void setID(uint32_t unitID, double timeSlice, int trackID, uint16_t depth=0)
Definition: CaloHitID.cc:41
uint32_t rotateUnitID(uint32_t, const G4Track *, const CastorShowerEvent &)
Definition: CastorSD.cc:299
CaloHitID currentID
Definition: CaloSD.h:142
G4LogicalVolume * lvCAST
Definition: CastorSD.h:49
CastorShowerLibrary * showerLibrary
Definition: CastorSD.h:47
virtual int setTrackID(const G4Step *)
Definition: CaloSD.cc:814
static bool isGammaElectronPositron(int pdgCode)
Geom::Theta< T > theta() const

◆ rotateUnitID()

uint32_t CastorSD::rotateUnitID ( uint32_t  unitID,
const G4Track *  track,
const CastorShowerEvent shower 
)
private

Definition at line 299 of file CastorSD.cc.

References printConversionInfo::aux, CastorShowerEvent::getPrimPhi(), createfilelist::int, M_PI, fileinputsource_cfi::sec, HLT_2022v12_cff::track, and reco::btau::trackPhi.

Referenced by getFromLibrary().

299  {
300  // ==============================================================
301  //
302  // o Exploit Castor phi symmetry to return newUnitID for
303  // shower hits based on track phi coordinate
304  //
305  // ==============================================================
306 
307  // Get 'track' phi:
308  double trackPhi = track->GetPosition().phi();
309  if (trackPhi < 0)
310  trackPhi += 2 * M_PI;
311  // Get phi from primary that gave rise to SL 'shower':
312  double showerPhi = shower.getPrimPhi();
313  if (showerPhi < 0)
314  showerPhi += 2 * M_PI;
315  // Delta phi:
316 
317  // Find the OctSector for which 'track' and 'shower' belong
318  int trackOctSector = (int)(trackPhi / (M_PI / 4));
319  int showerOctSector = (int)(showerPhi / (M_PI / 4));
320 
321  uint32_t newUnitID;
322  uint32_t sec = ((unitID >> 4) & 0xF);
323  uint32_t complement = (unitID & 0xFFFFFF0F);
324 
325  // Get 'track' z:
326  double trackZ = track->GetPosition().z();
327 
328  int aux;
329  int dSec = 2 * (trackOctSector - showerOctSector);
330  if (trackZ > 0) // Good for revision 1.9 of CastorNumberingScheme
331  {
332  int sec1 = sec - dSec;
333  // sec -= dSec ;
334  if (sec1 < 0)
335  sec1 += 16;
336  if (sec1 > 15)
337  sec1 -= 16;
338  sec = (uint32_t)(sec1);
339  } else {
340  if (dSec < 0)
341  sec += 16;
342  sec += dSec;
343  aux = (int)(sec / 16);
344  sec -= aux * 16;
345  }
346  sec = sec << 4;
347  newUnitID = complement | sec;
348 
349 #ifdef EDM_ML_DEBUG
350  if (newUnitID != unitID) {
351  edm::LogVerbatim("ForwardSim") << "\n CastorSD::rotateUnitID: "
352  << "\n unitID = " << unitID << "\n newUnitID = " << newUnitID;
353  }
354 #endif
355 
356  return newUnitID;
357 }
Log< level::Info, true > LogVerbatim
float getPrimPhi() const
#define M_PI

◆ setDetUnitId()

uint32_t CastorSD::setDetUnitId ( const G4Step *  step)
overridevirtual

Implements CaloSD.

Definition at line 283 of file CastorSD.cc.

References CastorNumberingScheme::getUnitID(), and numberingScheme.

283  {
284  return (numberingScheme == nullptr ? 0 : numberingScheme->getUnitID(aStep));
285 }
virtual uint32_t getUnitID(const G4Step *aStep) const
CastorNumberingScheme * numberingScheme
Definition: CastorSD.h:46

◆ setNumberingScheme()

void CastorSD::setNumberingScheme ( CastorNumberingScheme scheme)

Definition at line 289 of file CastorSD.cc.

References numberingScheme, and generator_cfi::scheme.

Referenced by CastorSD().

289  {
290  if (scheme != nullptr) {
291  edm::LogVerbatim("ForwardSim") << "CastorSD: updates numbering scheme for " << GetName();
292  delete numberingScheme;
294  }
295 }
Log< level::Info, true > LogVerbatim
CastorNumberingScheme * numberingScheme
Definition: CastorSD.h:46

Member Data Documentation

◆ energyThresholdSL

double CastorSD::energyThresholdSL
private

Definition at line 52 of file CastorSD.h.

Referenced by CastorSD(), and getFromLibrary().

◆ lvC3EF

G4LogicalVolume* CastorSD::lvC3EF
private

Definition at line 48 of file CastorSD.h.

Referenced by CastorSD(), and getEnergyDeposit().

◆ lvC3HF

G4LogicalVolume * CastorSD::lvC3HF
private

Definition at line 48 of file CastorSD.h.

Referenced by CastorSD(), and getEnergyDeposit().

◆ lvC4EF

G4LogicalVolume * CastorSD::lvC4EF
private

Definition at line 48 of file CastorSD.h.

Referenced by CastorSD(), and getEnergyDeposit().

◆ lvC4HF

G4LogicalVolume * CastorSD::lvC4HF
private

Definition at line 48 of file CastorSD.h.

Referenced by CastorSD(), and getEnergyDeposit().

◆ lvCAST

G4LogicalVolume* CastorSD::lvCAST
private

Definition at line 49 of file CastorSD.h.

Referenced by CastorSD(), and getFromLibrary().

◆ non_compensation_factor

double CastorSD::non_compensation_factor
private

Definition at line 53 of file CastorSD.h.

Referenced by CastorSD(), getEnergyDeposit(), and getFromLibrary().

◆ numberingScheme

CastorNumberingScheme* CastorSD::numberingScheme
private

Definition at line 46 of file CastorSD.h.

Referenced by setDetUnitId(), and setNumberingScheme().

◆ showerLibrary

CastorShowerLibrary* CastorSD::showerLibrary
private

Definition at line 47 of file CastorSD.h.

Referenced by CastorSD(), getFromLibrary(), and ~CastorSD().

◆ useShowerLibrary

bool CastorSD::useShowerLibrary
private

Definition at line 51 of file CastorSD.h.

Referenced by CastorSD(), and getFromLibrary().