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 DDCompactView &, 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 DDCompactView &cpv, 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
 
void PrintAll () override
 
G4bool ProcessHits (G4Step *step, G4TouchableHistory *) override
 
bool ProcessHits (G4GFlashSpot *aSpot, G4TouchableHistory *) override
 
 ~CaloSD () override
 
- Public Member Functions inherited from SensitiveCaloDetector
 SensitiveCaloDetector (const std::string &iname, const DDCompactView &cpv, const SensitiveDetectorCatalog &clg, edm::ParameterSet const &p)
 
- Public Member Functions inherited from SensitiveDetector
void EndOfEvent (G4HCofThisEvent *eventHC) override
 
const std::vector< std::string > & getNames () const
 
void Initialize (G4HCofThisEvent *eventHC) override
 
 SensitiveDetector (const std::string &iname, const DDCompactView &cpv, const SensitiveDetectorCatalog &, edm::ParameterSet const &p)
 
 ~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 *)
 
virtual void endEvent ()
 
virtual bool filterHit (CaloG4Hit *, double)
 
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 *)
 
virtual void initEvent (const BeginOfEvent *)
 
virtual void initRun ()
 
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)
 
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 }
 
- 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 30 of file CastorSD.h.

Constructor & Destructor Documentation

CastorSD::CastorSD ( const std::string &  name,
const DDCompactView cpv,
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(), GeV, lvC3EF, lvC3HF, lvC4EF, lvC4HF, lvCAST, non_compensation_factor, setNumberingScheme(), CaloSD::setParameterized(), showerLibrary, and useShowerLibrary.

34  :
35  CaloSD(name, cpv, clg, p, manager), numberingScheme(nullptr), lvC3EF(nullptr),
36  lvC3HF(nullptr), lvC4EF(nullptr), lvC4HF(nullptr), lvCAST(nullptr) {
37 
38  edm::ParameterSet m_CastorSD = p.getParameter<edm::ParameterSet>("CastorSD");
39  useShowerLibrary = m_CastorSD.getParameter<bool>("useShowerLibrary");
40  energyThresholdSL = m_CastorSD.getParameter<double>("minEnergyInGeVforUsingSLibrary");
41  energyThresholdSL = energyThresholdSL*GeV; // Convert GeV => MeV
42 
43  non_compensation_factor = m_CastorSD.getParameter<double>("nonCompensationFactor");
44 
45  if (useShowerLibrary) {
47  setParameterized(true);
48  }
50 
51  edm::LogInfo("ForwardSim")
52  << "********************************************************\n"
53  << "* Constructing a CastorSD with name " << GetName() << "\n"
54  << "* Using Castor Shower Library: " << useShowerLibrary << "\n"
55  << "********************************************************";
56 
57  const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
58  for (auto lv : *lvs) {
59  if (strcmp((lv->GetName()).c_str(),"C3EF") == 0) { lvC3EF = lv; }
60  if (strcmp((lv->GetName()).c_str(),"C3HF") == 0) { lvC3HF = lv; }
61  if (strcmp((lv->GetName()).c_str(),"C4EF") == 0) { lvC4EF = lv; }
62  if (strcmp((lv->GetName()).c_str(),"C4HF") == 0) { lvC4HF = lv; }
63  if (strcmp((lv->GetName()).c_str(),"CAST") == 0) { lvCAST = lv; }
64  if (lvC3EF != nullptr && lvC3HF != nullptr && lvC4EF != nullptr &&
65  lvC4HF != nullptr && lvCAST != nullptr) { break; }
66  }
67  edm::LogInfo("ForwardSim") << "CastorSD:: LogicalVolume pointers\n"
68  << lvC3EF << " for C3EF; " << lvC3HF
69  << " for C3HF; " << lvC4EF << " for C4EF; "
70  << lvC4HF << " for C4HF; "
71  << lvCAST << " for CAST. ";
72 }
T getParameter(std::string const &) const
const double GeV
Definition: MathUtil.h:16
double non_compensation_factor
Definition: CastorSD.h:55
CastorNumberingScheme * numberingScheme
Definition: CastorSD.h:48
bool useShowerLibrary
Definition: CastorSD.h:53
G4LogicalVolume * lvC4HF
Definition: CastorSD.h:50
G4LogicalVolume * lvC3HF
Definition: CastorSD.h:50
CaloSD(const std::string &aSDname, const DDCompactView &cpv, const SensitiveDetectorCatalog &clg, edm::ParameterSet const &p, const SimTrackManager *, float timeSlice=1., bool ignoreTkID=false)
Definition: CaloSD.cc:28
double energyThresholdSL
Definition: CastorSD.h:54
G4LogicalVolume * lvC4EF
Definition: CastorSD.h:50
G4LogicalVolume * lvCAST
Definition: CastorSD.h:51
CastorShowerLibrary * showerLibrary
Definition: CastorSD.h:49
void setNumberingScheme(CastorNumberingScheme *scheme)
Definition: CastorSD.cc:292
G4LogicalVolume * lvC3EF
Definition: CastorSD.h:50
void setParameterized(bool val)
Definition: CaloSD.h:99
CastorSD::~CastorSD ( )
override

Definition at line 76 of file CastorSD.cc.

References showerLibrary.

76  {
77  delete showerLibrary;
78 }
CastorShowerLibrary * showerLibrary
Definition: CastorSD.h:49

Member Function Documentation

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

Reimplemented from CaloSD.

Definition at line 82 of file CastorSD.cc.

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

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

Reimplemented from CaloSD.

Definition at line 363 of file CastorSD.cc.

References funct::abs(), pfBoostedDoubleSVAK8TagInfos_cfi::beta, ALCARECOTkAlJpsiMuMu_cff::charge, CaloSD::currentID, dot(), CaloSD::edepositEM, CaloSD::edepositHAD, energyThresholdSL, CastorShowerEvent::getDetID(), CastorShowerEvent::getNhit(), CastorShowerEvent::getNphotons(), CastorShowerEvent::getPrimE(), CastorShowerLibrary::getShowerHits(), CastorShowerEvent::getTime(), GeV, TrackInformation::hasCastorHit(), hfClusterShapes_cfi::hits, mps_fire::i, G4TrackToParticleID::isGammaElectronPositron(), G4TrackToParticleID::isStableHadronIon(), lvC3EF, lvC3HF, lvC4EF, lvC4HF, lvCAST, M_PI, non_compensation_factor, CaloSD::processHit(), dttmaxenums::R, CaloSD::resetForNewPrimary(), rotateUnitID(), Scenarios_cff::scale, TrackInformation::setCastorHitPID(), CaloHitID::setID(), CaloSD::setTrackID(), showerLibrary, mathSSE::sqrt(), theta(), ntuplemaker::time, and useShowerLibrary.

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

Definition at line 304 of file CastorSD.cc.

References printConversionInfo::aux, CastorShowerEvent::getPrimPhi(), createfilelist::int, LogDebug, M_PI, and reco::btau::trackPhi.

Referenced by getFromLibrary().

304  {
305 // ==============================================================
306 //
307 // o Exploit Castor phi symmetry to return newUnitID for
308 // shower hits based on track phi coordinate
309 //
310 // ==============================================================
311 
312  // Get 'track' phi:
313  double trackPhi = track->GetPosition().phi();
314  if(trackPhi<0) trackPhi += 2*M_PI ;
315  // Get phi from primary that gave rise to SL 'shower':
316  double showerPhi = shower.getPrimPhi();
317  if(showerPhi<0) showerPhi += 2*M_PI ;
318  // Delta phi:
319 
320  // Find the OctSector for which 'track' and 'shower' belong
321  int trackOctSector = (int) ( trackPhi / (M_PI/4) ) ;
322  int showerOctSector = (int) ( showerPhi / (M_PI/4) ) ;
323 
324  uint32_t newUnitID;
325  uint32_t sec = ( ( unitID>>4 ) & 0xF ) ;
326  uint32_t complement = ( unitID & 0xFFFFFF0F ) ;
327 
328  // Get 'track' z:
329  double trackZ = track->GetPosition().z();
330 
331  int aux ;
332  int dSec = 2*(trackOctSector - showerOctSector) ;
333  if(trackZ>0) // Good for revision 1.9 of CastorNumberingScheme
334  {
335  int sec1 = sec-dSec;
336  // sec -= dSec ;
337  if(sec1<0) sec1 += 16;
338  if(sec1>15) sec1 -= 16;
339  sec = (uint32_t)(sec1);
340  } else {
341  if( dSec<0 ) 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 debugLog
350  if(newUnitID != unitID) {
351  LogDebug("ForwardSim") << "\n CastorSD::rotateUnitID: "
352  << "\n unitID = " << unitID
353  << "\n newUnitID = " << newUnitID ;
354  }
355 #endif
356 
357  return newUnitID ;
358 
359 }
#define LogDebug(id)
#define M_PI
float getPrimPhi() const
uint32_t CastorSD::setDetUnitId ( const G4Step *  step)
overridevirtual

Implements CaloSD.

Definition at line 286 of file CastorSD.cc.

References CastorNumberingScheme::getUnitID(), and numberingScheme.

286  {
287  return (numberingScheme == nullptr ? 0 : numberingScheme->getUnitID(aStep));
288 }
CastorNumberingScheme * numberingScheme
Definition: CastorSD.h:48
virtual uint32_t getUnitID(const G4Step *aStep) const
void CastorSD::setNumberingScheme ( CastorNumberingScheme scheme)

Definition at line 292 of file CastorSD.cc.

References numberingScheme.

Referenced by CastorSD().

292  {
293 
294  if (scheme != nullptr) {
295  edm::LogInfo("ForwardSim") << "CastorSD: updates numbering scheme for "
296  << GetName();
297  delete numberingScheme;
298  numberingScheme = scheme;
299  }
300 }
CastorNumberingScheme * numberingScheme
Definition: CastorSD.h:48

Member Data Documentation

double CastorSD::energyThresholdSL
private

Definition at line 54 of file CastorSD.h.

Referenced by CastorSD(), and getFromLibrary().

G4LogicalVolume* CastorSD::lvC3EF
private

Definition at line 50 of file CastorSD.h.

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

G4LogicalVolume * CastorSD::lvC3HF
private

Definition at line 50 of file CastorSD.h.

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

G4LogicalVolume * CastorSD::lvC4EF
private

Definition at line 50 of file CastorSD.h.

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

G4LogicalVolume * CastorSD::lvC4HF
private

Definition at line 50 of file CastorSD.h.

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

G4LogicalVolume* CastorSD::lvCAST
private

Definition at line 51 of file CastorSD.h.

Referenced by CastorSD(), and getFromLibrary().

double CastorSD::non_compensation_factor
private

Definition at line 55 of file CastorSD.h.

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

CastorNumberingScheme* CastorSD::numberingScheme
private

Definition at line 48 of file CastorSD.h.

Referenced by setDetUnitId(), and setNumberingScheme().

CastorShowerLibrary* CastorSD::showerLibrary
private

Definition at line 49 of file CastorSD.h.

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

bool CastorSD::useShowerLibrary
private

Definition at line 53 of file CastorSD.h.

Referenced by CastorSD(), and getFromLibrary().