CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 (G4String, const DDCompactView &, SensitiveDetectorCatalog &clg, edm::ParameterSet const &, const SimTrackManager *)
 
virtual double getEnergyDeposit (G4Step *)
 
virtual uint32_t setDetUnitId (G4Step *step)
 
void setNumberingScheme (CastorNumberingScheme *scheme)
 
virtual ~CastorSD ()
 
- Public Member Functions inherited from CaloSD
 CaloSD (G4String aSDname, const DDCompactView &cpv, SensitiveDetectorCatalog &clg, edm::ParameterSet const &p, const SimTrackManager *, int tSlice=1, bool ignoreTkID=false)
 
virtual void clear ()
 
virtual void DrawAll ()
 
virtual void EndOfEvent (G4HCofThisEvent *eventHC)
 
void fillHits (edm::PCaloHitContainer &, std::string n)
 
virtual void Initialize (G4HCofThisEvent *HCE)
 
virtual void PrintAll ()
 
virtual bool ProcessHits (G4Step *step, G4TouchableHistory *tHistory)
 
virtual bool ProcessHits (G4GFlashSpot *aSpot, G4TouchableHistory *)
 
virtual ~CaloSD ()
 
- Public Member Functions inherited from SensitiveCaloDetector
 SensitiveCaloDetector (std::string &iname, const DDCompactView &cpv, SensitiveDetectorCatalog &clg, edm::ParameterSet const &p)
 
- Public Member Functions inherited from SensitiveDetector
virtual void AssignSD (std::string &vname)
 
Local3DPoint ConvertToLocal3DPoint (const G4ThreeVector &point)
 
Local3DPoint FinalStepPosition (G4Step *s, coordinates)
 
virtual std::vector< std::string > getNames ()
 
Local3DPoint InitialStepPosition (G4Step *s, coordinates)
 
std::string nameOfSD ()
 
void NaNTrap (G4Step *step)
 
void Register ()
 
 SensitiveDetector (std::string &iname, const DDCompactView &cpv, SensitiveDetectorCatalog &, edm::ParameterSet const &p)
 
virtual ~SensitiveDetector ()
 
- 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

virtual void initRun ()
 
- Protected Member Functions inherited from CaloSD
G4bool checkHit ()
 
virtual void clearHits ()
 
CaloG4HitcreateNewHit ()
 
virtual bool filterHit (CaloG4Hit *, double)
 
double getAttenuation (G4Step *aStep, double birk1, double birk2, double birk3)
 
virtual uint16_t getDepth (G4Step *)
 
int getNumberOfHits ()
 
double getResponseWt (G4Track *)
 
virtual G4bool getStepInfo (G4Step *aStep)
 
virtual int getTrackID (G4Track *)
 
G4bool hitExists ()
 
void resetForNewPrimary (const G4ThreeVector &, double)
 
G4ThreeVector setToGlobal (const G4ThreeVector &, const G4VTouchable *)
 
G4ThreeVector setToLocal (const G4ThreeVector &, const G4VTouchable *)
 
virtual void update (const BeginOfRun *)
 This routine will be called when the appropriate signal arrives. More...
 
virtual void update (const BeginOfEvent *)
 This routine will be called when the appropriate signal arrives. More...
 
virtual void update (const BeginOfTrack *trk)
 This routine will be called when the appropriate signal arrives. More...
 
virtual void update (const EndOfTrack *trk)
 This routine will be called when the appropriate signal arrives. More...
 
virtual void update (const ::EndOfEvent *)
 
void updateHit (CaloG4Hit *)
 
- 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

void getFromLibrary (G4Step *)
 
uint32_t rotateUnitID (uint32_t, G4Track *, const CastorShowerEvent &)
 
int setTrackID (G4Step *)
 

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

- Public Types inherited from SensitiveDetector
enum  coordinates { WorldCoordinates, LocalCoordinates }
 
- Protected Attributes inherited from CaloSD
int checkHits
 
double correctT
 
bool corrTOFBeam
 
CaloG4HitcurrentHit
 
CaloHitID currentID
 
float edepositEM
 
float edepositHAD
 
double eminHit
 
double eminHitD
 
G4int emPDG
 
double energyCut
 
G4ThreeVector entranceLocal
 
G4ThreeVector entrancePoint
 
G4int epPDG
 
bool forceSave
 
G4int gammaPDG
 
float incidentEnergy
 
double kmaxIon
 
double kmaxNeutron
 
double kmaxProton
 
const SimTrackManagerm_trackManager
 
G4ThreeVector posGlobal
 
G4StepPoint * preStepPoint
 
CaloHitID previousID
 
int primIDSaved
 
bool runInit
 
bool suppressHeavy
 
G4Track * theTrack
 
double tmaxHit
 
bool useMap
 

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 ( G4String  name,
const DDCompactView cpv,
SensitiveDetectorCatalog clg,
edm::ParameterSet const &  p,
const SimTrackManager manager 
)

Definition at line 31 of file CastorSD.cc.

References energyThresholdSL, edm::ParameterSet::getParameter(), lvC3EF, lvC3HF, lvC4EF, lvC4HF, lvCAST, non_compensation_factor, setNumberingScheme(), showerLibrary, and useShowerLibrary.

34  :
35  CaloSD(name, cpv, clg, p, manager), numberingScheme(0), lvC3EF(0),
36  lvC3HF(0), lvC4EF(0), lvC4HF(0), lvCAST(0) {
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) showerLibrary = new CastorShowerLibrary(name, p);
46 
48 
49  edm::LogInfo("ForwardSim")
50  << "***************************************************\n"
51  << "* *\n"
52  << "* Constructing a CastorSD with name " << GetName() << "\n"
53  << "* *\n"
54  << "***************************************************";
55 
56  const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
57  std::vector<G4LogicalVolume*>::const_iterator lvcite;
58  for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
59  if (strcmp(((*lvcite)->GetName()).c_str(),"C3EF") == 0) lvC3EF = (*lvcite);
60  if (strcmp(((*lvcite)->GetName()).c_str(),"C3HF") == 0) lvC3HF = (*lvcite);
61  if (strcmp(((*lvcite)->GetName()).c_str(),"C4EF") == 0) lvC4EF = (*lvcite);
62  if (strcmp(((*lvcite)->GetName()).c_str(),"C4HF") == 0) lvC4HF = (*lvcite);
63  if (strcmp(((*lvcite)->GetName()).c_str(),"CAST") == 0) lvCAST = (*lvcite);
64  if (lvC3EF != 0 && lvC3HF != 0 && lvC4EF != 0 && lvC4HF != 0 && lvCAST != 0) break;
65  }
66  edm::LogInfo("ForwardSim") << "CastorSD:: LogicalVolume pointers\n"
67  << lvC3EF << " for C3EF; " << lvC3HF
68  << " for C3HF; " << lvC4EF << " for C4EF; "
69  << lvC4HF << " for C4HF; "
70  << lvCAST << " for CAST. " << std::endl;
71 
72  // if(useShowerLibrary) edm::LogInfo("ForwardSim") << "\n Using Castor Shower Library \n";
73 
74 }
T getParameter(std::string const &) const
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
CaloSD(G4String aSDname, const DDCompactView &cpv, SensitiveDetectorCatalog &clg, edm::ParameterSet const &p, const SimTrackManager *, int tSlice=1, bool ignoreTkID=false)
Definition: CaloSD.cc:24
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:515
G4LogicalVolume * lvC3EF
Definition: CastorSD.h:48
CastorSD::~CastorSD ( )
virtual

Definition at line 78 of file CastorSD.cc.

References showerLibrary, and useShowerLibrary.

78  {
79  if (useShowerLibrary) delete showerLibrary;
80 }
bool useShowerLibrary
Definition: CastorSD.h:51
CastorShowerLibrary * showerLibrary
Definition: CastorSD.h:47

Member Function Documentation

double CastorSD::getEnergyDeposit ( G4Step *  aStep)
virtual

Reimplemented from CaloSD.

Definition at line 95 of file CastorSD.cc.

References a, funct::abs(), beta, DeDxDiscriminatorTools::charge(), dot(), energyThresholdSL, eta(), getFromLibrary(), create_public_lumi_plots::log, LogDebug, lvC3EF, lvC3HF, lvC4EF, lvC4HF, lvCAST, M_PI, max(), bookConverter::min, SensitiveDetector::name, NULL, phi, pi, CaloSD::preStepPoint, dttmaxenums::R, alignCSCRings::r, mathSSE::sqrt(), AlCaHLTBitMon_QueryRunRegistry::string, funct::tan(), theta(), CaloSD::theTrack, and useShowerLibrary.

95  {
96 
97  float NCherPhot = 0.;
98 
99  if (aStep == NULL) {
100  return 0;
101  } else {
102  // preStepPoint information *********************************************
103 
104  G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
105  G4VPhysicalVolume* currentPV = preStepPoint->GetPhysicalVolume();
106  G4LogicalVolume* currentLV = currentPV->GetLogicalVolume();
107  G4String name = currentPV->GetName();
108  std::string nameVolume;
109  nameVolume.assign(name,0,4);
110 
111 #ifdef debugLog
112  G4SteppingControl stepControlFlag = aStep->GetControlFlag();
113  if (aStep->IsFirstStepInVolume())
114  LogDebug("ForwardSim") << "CastorSD::getEnergyDeposit:"
115  << "\n IsFirstStepInVolume " ;
116 #endif
117 
118  // Get theTrack
119  G4Track* theTrack = aStep->GetTrack();
120 
121 #ifdef debugLog
122  if (useShowerLibrary && currentLV==lvCAST) {
123  LogDebug("ForwardSim") << "CastorSD::getEnergyDeposit:"
124  << "\n TrackID , ParentID , ParticleName ,"
125  << " eta , phi , z , time ,"
126  << " K , E , Mom "
127  << "\n TRACKINFO: "
128  << theTrack->GetTrackID()
129  << " , "
130  << theTrack->GetParentID()
131  << " , "
132  << theTrack->GetDefinition()->GetParticleName()
133  << " , "
134  << theTrack->GetPosition().eta()
135  << " , "
136  << theTrack->GetPosition().phi()
137  << " , "
138  << theTrack->GetPosition().z()
139  << " , "
140  << theTrack->GetGlobalTime()
141  << " , "
142  << theTrack->GetKineticEnergy()
143  << " , "
144  << theTrack->GetTotalEnergy()
145  << " , "
146  << theTrack->GetMomentum().mag() ;
147  if(theTrack->GetTrackID() != 1)
148  LogDebug("ForwardSim") << "CastorSD::getEnergyDeposit:"
149  << "\n CurrentStepNumber , TrackID , Particle , VertexPosition ,"
150  << " LogicalVolumeAtVertex , CreatorProcess"
151  << "\n TRACKINFO2: "
152  << theTrack->GetCurrentStepNumber()
153  << " , "
154  << theTrack->GetTrackID()
155  << " , "
156  << theTrack->GetDefinition()->GetParticleName()
157  << " , "
158  << theTrack->GetVertexPosition()
159  << " , "
160  << theTrack->GetLogicalVolumeAtVertex()->GetName()
161  << " , "
162  << theTrack->GetCreatorProcess()->GetProcessName() ;
163  } // end of if(useShowerLibrary)
164 #endif
165 
166  // if particle moves from interaction point or "backwards (halo)
167  bool backward = false;
168  G4ThreeVector hitPoint = preStepPoint->GetPosition();
169  G4ThreeVector hit_mom = preStepPoint->GetMomentumDirection();
170  double zint = hitPoint.z();
171  double pz = hit_mom.z();
172 
173  // Check if theTrack moves backward
174  if (pz * zint < 0.) backward = true;
175 
176  // Check that theTrack is above the energy threshold to use Shower Library
177  bool aboveThreshold = false;
178  if(theTrack->GetKineticEnergy() > energyThresholdSL) aboveThreshold = true;
179 
180  // Check if theTrack is a muon (if so, DO NOT use Shower Library)
181  bool notaMuon = true;
182  G4int mumPDG = 13;
183  G4int mupPDG = -13;
184  G4int parCode = theTrack->GetDefinition()->GetPDGEncoding();
185  if (parCode == mupPDG || parCode == mumPDG ) notaMuon = false;
186 
187  // angle condition
188  double theta_max = M_PI - 3.1305; // angle in radians corresponding to -5.2 eta
189  double R_mom=sqrt(hit_mom.x()*hit_mom.x() + hit_mom.y()*hit_mom.y());
190  double theta = atan2(R_mom,std::abs(pz));
191  bool angleok = false;
192  if ( theta < theta_max) angleok = true;
193 
194  // OkToUse
195  double R = sqrt(hitPoint.x()*hitPoint.x() + hitPoint.y()*hitPoint.y());
196  bool dot = false;
197  if ( zint < -14450. && R < 45.) dot = true;
198  bool inRange = true;
199  if ( zint < -14700. || R > 193.) inRange = false;
200  bool OkToUse = false;
201  if ( inRange && !dot) OkToUse = true;
202 
203  if (useShowerLibrary && aboveThreshold && notaMuon && (!backward) && OkToUse && angleok && currentLV == lvCAST ) {
204  // Use Castor shower library if energy is above threshold, is not a muon
205  // and is not moving backward
206  getFromLibrary(aStep);
207 
208 #ifdef debugLog
209  LogDebug("ForwardSim") << " Current logical volume is " << nameVolume ;
210 #endif
211  return NCherPhot;
212  } else {
213  // Usual calculations
214  // G4ThreeVector hitPoint = preStepPoint->GetPosition();
215  // G4ThreeVector hit_mom = preStepPoint->GetMomentumDirection();
216  G4double stepl = aStep->GetStepLength()/cm;
217  G4double beta = preStepPoint->GetBeta();
218  G4double charge = preStepPoint->GetCharge();
219  // G4VProcess* curprocess = preStepPoint->GetProcessDefinedStep();
220  // G4String namePr = preStepPoint->GetProcessDefinedStep()->GetProcessName();
221  // std::string nameProcess;
222  // nameProcess.assign(namePr,0,4);
223 
224  // G4LogicalVolume* lv = currentPV->GetLogicalVolume();
225  // G4Material* mat = lv->GetMaterial();
226  // G4double rad = mat->GetRadlen();
227 
228 
229  // postStepPoint information *********************************************
230  G4StepPoint* postStepPoint= aStep->GetPostStepPoint();
231  G4VPhysicalVolume* postPV = postStepPoint->GetPhysicalVolume();
232  G4String postname = postPV->GetName();
233  std::string postnameVolume;
234  postnameVolume.assign(postname,0,4);
235 
236  // theTrack information *************************************************
237  // G4Track* theTrack = aStep->GetTrack();
238  //G4double entot = theTrack->GetTotalEnergy();
239  G4ThreeVector vert_mom = theTrack->GetVertexMomentumDirection();
240 
241  G4ThreeVector localPoint = theTrack->GetTouchable()->GetHistory()->
242  GetTopTransform().TransformPoint(hitPoint);
243 
244  // calculations... *************************************************
245  float phi = -100.;
246  if (vert_mom.x() != 0) phi = atan2(vert_mom.y(),vert_mom.x());
247  if (phi < 0.) phi += twopi;
248  G4String particleType = theTrack->GetDefinition()->GetParticleName();
249 #ifdef debugLog
250  float costheta =vert_mom.z()/sqrt(vert_mom.x()*vert_mom.x()+
251  vert_mom.y()*vert_mom.y()+
252  vert_mom.z()*vert_mom.z());
253  float theta = acos(std::min(std::max(costheta,float(-1.)),float(1.)));
254  float eta = -log(tan(theta/2));
255  G4int primaryID = theTrack->GetTrackID();
256  // *************************************************
257 
258 
259  // *************************************************
260  double edep = aStep->GetTotalEnergyDeposit();
261 #endif
262 
263  // *************************************************
264 
265 
266  // *************************************************
267  // take into account light collection curve for plate
268  // double weight = curve_Castor(nameVolume, preStepPoint);
269  // double edep = aStep->GetTotalEnergyDeposit() * weight;
270  // *************************************************
271 
272 
273  // *************************************************
274  /* comments for sensitive volumes:
275  C001 ...-... CP06,CBU1 ...-...CALM --- > fibres and bundle
276  for first release of CASTOR
277  CASF --- > quartz plate for first and second releases of CASTOR
278  GF2Q, GFNQ, GR2Q, GRNQ
279  for tests with my own test geometry of HF (on ask of Gavrilov)
280  C3TF, C4TF - for third release of CASTOR
281  */
282  double meanNCherPhot=0;
283 
284  if (currentLV == lvC3EF || currentLV == lvC4EF || currentLV == lvC3HF ||
285  currentLV == lvC4HF) {
286  // if(nameVolume == "C3EF" || nameVolume == "C4EF" || nameVolume == "C3HF" || nameVolume == "C4HF") {
287 
288  float bThreshold = 0.67;
289  float nMedium = 1.4925;
290  // float photEnSpectrDL = (1./400.nm-1./700.nm)*10000000.cm/nm; /* cm-1 */
291  // float photEnSpectrDL = 10714.285714;
292 
293  float photEnSpectrDE = 1.24; /* see below */
294  /* E = 2pi*(1./137.)*(eV*cm/370.)/lambda = */
295  /* = 12.389184*(eV*cm)/lambda */
296  /* Emax = 12.389184*(eV*cm)/400nm*10-7cm/nm = 3.01 eV */
297  /* Emin = 12.389184*(eV*cm)/700nm*10-7cm/nm = 1.77 eV */
298  /* delE = Emax - Emin = 1.24 eV --> */
299  /* */
300  /* default for Castor nameVolume == "CASF" or (C3TF & C4TF) */
301 
302  float thFullRefl = 23.; /* 23.dergee */
303  float thFullReflRad = thFullRefl*pi/180.;
304 
305  /* default for Castor nameVolume == "CASF" or (C3TF & C4TF) */
306  float thFibDir = 45.; /* .dergee */
307  /* for test HF geometry volumes:
308  if(nameVolume == "GF2Q" || nameVolume == "GFNQ" ||
309  nameVolume == "GR2Q" || nameVolume == "GRNQ")
310  thFibDir = 0.0; // .dergee
311  */
312  float thFibDirRad = thFibDir*pi/180.;
313  /* */
314  /* */
315 
316  // at which theta the point is located:
317  // float th1 = hitPoint.theta();
318 
319  // theta of charged particle in LabRF(hit momentum direction):
320  float costh =hit_mom.z()/sqrt(hit_mom.x()*hit_mom.x()+
321  hit_mom.y()*hit_mom.y()+
322  hit_mom.z()*hit_mom.z());
323  if (zint < 0) costh = -costh;
324  float th = acos(std::min(std::max(costh,float(-1.)),float(1.)));
325 
326  // just in case (can do bot use):
327  if (th < 0.) th += twopi;
328 
329 
330 
331  // theta of cone with Cherenkov photons w.r.t.direction of charged part.:
332  float costhcher =1./(nMedium*beta);
333  float thcher = acos(std::min(std::max(costhcher,float(-1.)),float(1.)));
334 
335  // diff thetas of charged part. and quartz direction in LabRF:
336  float DelFibPart = fabs(th - thFibDirRad);
337 
338  // define real distances:
339  float d = fabs(tan(th)-tan(thFibDirRad));
340 
341  // float a = fabs(tan(thFibDirRad)-tan(thFibDirRad+thFullReflRad));
342  // float r = fabs(tan(th)-tan(th+thcher));
343 
344  float a = tan(thFibDirRad)+tan(fabs(thFibDirRad-thFullReflRad));
345  float r = tan(th)+tan(fabs(th-thcher));
346 
347 
348  // define losses d_qz in cone of full reflection inside quartz direction
349  float d_qz;
350 #ifdef debugLog
351  float variant;
352 #endif
353  if(DelFibPart > (thFullReflRad + thcher) ) {
354  d_qz = 0.;
355 #ifdef debugLog
356  variant=0.;
357 #endif
358  }
359  // if(d > (r+a) ) {d_qz = 0.; variant=0.;}
360  else {
361  if((th + thcher) < (thFibDirRad+thFullReflRad) &&
362  (th - thcher) > (thFibDirRad-thFullReflRad)) {
363  d_qz = 1.;
364 #ifdef debugLog
365  variant=1.;
366 #endif
367  }
368  // if((DelFibPart + thcher) < thFullReflRad ) {d_qz = 1.; variant=1.;}
369  // if((d+r) < a ) {d_qz = 1.; variant=1.;}
370  else {
371  if((thFibDirRad + thFullReflRad) < (th + thcher) &&
372  (thFibDirRad - thFullReflRad) > (th - thcher) ) {
373  // if((thcher - DelFibPart ) > thFullReflRad )
374  // if((r-d ) > a )
375  d_qz = 0.;
376 #ifdef debugLog
377  variant=2.;
378 #endif
379  } else {
380  // if((thcher + DelFibPart ) > thFullReflRad &&
381  // thcher < (DelFibPart+thFullReflRad) )
382  // {
383  d_qz = 0.;
384 #ifdef debugLog
385  variant=3.;
386 #endif
387 
388  // use crossed length of circles(cone projection)
389  // dC1/dC2 :
390  float arg_arcos = 0.;
391  float tan_arcos = 2.*a*d;
392  if(tan_arcos != 0.) arg_arcos =(r*r-a*a-d*d)/tan_arcos;
393  arg_arcos = fabs(arg_arcos);
394  float th_arcos = acos(std::min(std::max(arg_arcos,float(-1.)),float(1.)));
395  d_qz = th_arcos/pi/2.;
396  d_qz = fabs(d_qz);
397 
398 
399 
400  // }
401  // else
402  // {
403  // d_qz = 0.; variant=4.;
404  //#ifdef debugLog
405  // std::cout <<" ===============>variant 4 information: <===== " <<std::endl;
406  // std::cout <<" !!!!!!!!!!!!!!!!!!!!!! variant = " << variant <<std::endl;
407  //#endif
408  //
409  // }
410  }
411  }
412  }
413 
414 
415  // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
416 
417  if(charge != 0. && beta > bThreshold ) {
418 
419  meanNCherPhot = 370.*charge*charge*
420  ( 1. - 1./(nMedium*nMedium*beta*beta) )*
421  photEnSpectrDE*stepl;
422 
423  G4int poissNCherPhot = (G4int) G4Poisson(meanNCherPhot);
424 
425  if(poissNCherPhot < 0) poissNCherPhot = 0;
426 
427  float effPMTandTransport = 0.19;
428  double ReflPower = 0.1;
429  double proba = d_qz + (1-d_qz)*ReflPower;
430  NCherPhot = poissNCherPhot*effPMTandTransport*proba*0.307;
431 
432 
433 #ifdef debugLog
434  float thgrad = th*180./pi;
435  float thchergrad = thcher*180./pi;
436  float DelFibPartgrad = DelFibPart*180./pi;
437  LogDebug("ForwardSim") << " ==============================> start all "
438  << "information:<========= \n" << " =====> for "
439  << "test:<=== \n" << " variant = " << variant
440  << "\n thgrad = " << thgrad << "\n thchergrad "
441  << "= " << thchergrad << "\n DelFibPartgrad = "
442  << DelFibPartgrad << "\n d_qz = " << d_qz
443  << "\n =====> Start Step Information <=== \n"
444  << " ===> calo preStepPoint info <=== \n"
445  << " hitPoint = " << hitPoint << "\n"
446  << " hitMom = " << hit_mom << "\n"
447  << " stepControlFlag = " << stepControlFlag
448  // << "\n curprocess = " << curprocess << "\n"
449  // << " nameProcess = " << nameProcess
450  << "\n charge = " << charge << "\n"
451  << " beta = " << beta << "\n"
452  << " bThreshold = " << bThreshold << "\n"
453  << " thgrad =" << thgrad << "\n"
454  << " effPMTandTransport=" << effPMTandTransport
455  // << "\n volume = " << name
456  << "\n nameVolume = " << nameVolume << "\n"
457  << " nMedium = " << nMedium << "\n"
458  // << " rad length = " << rad << "\n"
459  // << " material = " << mat << "\n"
460  << " stepl = " << stepl << "\n"
461  << " photEnSpectrDE = " << photEnSpectrDE <<"\n"
462  << " edep = " << edep << "\n"
463  << " ===> calo theTrack info <=== " << "\n"
464  << " particleType = " << particleType << "\n"
465  << " primaryID = " << primaryID << "\n"
466  << " entot= " << theTrack->GetTotalEnergy() << "\n"
467  << " vert_eta= " << eta << "\n"
468  << " vert_phi= " << phi << "\n"
469  << " vert_mom= " << vert_mom << "\n"
470  << " ===> calo hit preStepPointinfo <=== "<<"\n"
471  << " local point = " << localPoint << "\n"
472  << " ==============================> final info"
473  << ": <=== \n"
474  << " meanNCherPhot = " << meanNCherPhot << "\n"
475  << " poissNCherPhot = " << poissNCherPhot <<"\n"
476  << " NCherPhot = " << NCherPhot;
477 #endif
478 
479  // Included by WC
480  // std::cout << "\n volume = " << name
481  // << "\n nameVolume = " << nameVolume << "\n"
482  // << "\n postvolume = " << postname
483  // << "\n postnameVolume = " << postnameVolume << "\n"
484  // << "\n particleType = " << particleType
485  // << "\n primaryID = " << primaryID << "\n";
486 
487  }
488  }
489 
490 
491 #ifdef debugLog
492  LogDebug("ForwardSim") << "CastorSD:: " << nameVolume
493  // << " Light Collection Efficiency " << weight
494  << " Weighted Energy Deposit " << edep/MeV
495  << " MeV\n";
496 #endif
497  // Temporary member for testing purpose only...
498  // unit_id = setDetUnitId(aStep);
499  // if(NCherPhot != 0) std::cout << "\n UnitID = " << unit_id << " ; NCherPhot = " << NCherPhot ;
500 
501  return NCherPhot;
502  }
503  }
504  return 0;
505 }
#define LogDebug(id)
const double beta
Geom::Theta< T > theta() const
#define NULL
Definition: scimark2.h:8
bool useShowerLibrary
Definition: CastorSD.h:51
T eta() const
double charge(const std::vector< uint8_t > &Ampls)
void getFromLibrary(G4Step *)
Definition: CastorSD.cc:615
G4LogicalVolume * lvC4HF
Definition: CastorSD.h:48
const T & max(const T &a, const T &b)
T sqrt(T t)
Definition: SSEVec.h:48
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
G4Track * theTrack
Definition: CaloSD.h:118
double energyThresholdSL
Definition: CastorSD.h:52
G4StepPoint * preStepPoint
Definition: CaloSD.h:120
G4LogicalVolume * lvC4EF
Definition: CastorSD.h:48
#define M_PI
Definition: BFit3D.cc:3
G4LogicalVolume * lvCAST
Definition: CastorSD.h:49
G4LogicalVolume * lvC3EF
Definition: CastorSD.h:48
T dot(const Basic3DVector &v) const
Scalar product, or &quot;dot&quot; product, with a vector of same type.
double a
Definition: hdecay.h:121
double pi
Definition: DDAxes.h:10
void CastorSD::getFromLibrary ( G4Step *  aStep)
private

Definition at line 615 of file CastorSD.cc.

References CaloSD::checkHit(), CaloSD::createNewHit(), CaloSD::currentHit, CaloSD::currentID, CaloSD::edepositEM, CaloSD::edepositHAD, CaloSD::emPDG, CaloSD::epPDG, CaloSD::gammaPDG, CastorShowerEvent::getDetID(), CastorShowerEvent::getNhit(), CastorShowerEvent::getNphotons(), CastorShowerEvent::getPrimE(), CastorShowerLibrary::getShowerHits(), CastorShowerEvent::getTime(), GetVolume(), i, GetRecoTauVFromDQM_MC_cff::kk, LogDebug, non_compensation_factor, convertSQLiteXML::ok, CaloSD::posGlobal, CaloSD::preStepPoint, CaloSD::previousID, CaloSD::resetForNewPrimary(), rotateUnitID(), pileupReCalc_HLTpaths::scale, CaloHitID::setID(), setTrackID(), showerLibrary, CaloSD::theTrack, cond::rpcobgas::time, funct::true, and CaloSD::updateHit().

Referenced by getEnergyDeposit().

615  {
616 
618 //
619 // Method to get hits from the Shower Library
620 //
621 // CastorShowerEvent hits returned by getShowerHits are used to
622 // replace the full simulation of the shower from theTrack
623 //
624 // "updateHit" save the Hits to a CaloG4Hit container
625 //
627 
628  preStepPoint = aStep->GetPreStepPoint();
629  theTrack = aStep->GetTrack();
630  bool ok;
631 
632  // **** Call method to retrieve hits from the ShowerLibrary ****
633  CastorShowerEvent hits = showerLibrary->getShowerHits(aStep, ok);
634 
635  double etrack = preStepPoint->GetKineticEnergy();
636  int primaryID = setTrackID(aStep);
637  // int primaryID = theTrack->GetTrackID();
638 
639  // Reset entry point for new primary
640  posGlobal = preStepPoint->GetPosition();
641  resetForNewPrimary(posGlobal, etrack);
642 
643  // Check whether track is EM or HAD
644  G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
645  bool isEM , isHAD ;
646  if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG) {
647  isEM = true ; isHAD = false;
648  } else {
649  isEM = false; isHAD = true ;
650  }
651 
652 #ifdef debugLog
653  // edm::LogInfo("ForwardSim") << "\n CastorSD::getFromLibrary: "
654  LogDebug("ForwardSim") << "\n CastorSD::getFromLibrary: "
655  << hits.getNhit() << " hits for " << GetName() << " from "
656  << theTrack->GetDefinition()->GetParticleName() << " of "
657  << preStepPoint->GetKineticEnergy()/GeV << " GeV and trackID "
658  << theTrack->GetTrackID() ;
659 #endif
660 
661  // Scale to correct energy
662  double E_track = preStepPoint->GetTotalEnergy() ;
663  double E_SLhit = hits.getPrimE() * GeV ;
664  double scale = E_track/E_SLhit ;
665 
666  //Non compensation
667  if (isHAD){
668  scale=scale*non_compensation_factor; // if hadronic extend the scale with the non-compensation factor
669  } else {
670  scale=scale; // if electromagnetic, don't do anything
671  }
672 
673 
674 /* double theTrackEnergy = theTrack->GetTotalEnergy() ;
675 
676  if(fabs(theTrackEnergy-E_track)>10.) {
677  edm::LogInfo("ForwardSim") << "\n TrackID = " << theTrack->GetTrackID()
678  << "\n theTrackEnergy = " << theTrackEnergy
679  << "\n preStepPointEnergy = " << E_track ;
680  G4TrackVector tsec = *(aStep->GetSecondary());
681  for (unsigned int kk=0; kk<tsec.size(); kk++) {
682  edm::LogInfo("ForwardSim") << "CastorSD::getFromLibrary:"
683  << "\n tsec[" << kk << "]->GetTrackID() = "
684  << tsec[kk]->GetTrackID()
685  << " with energy "
686  << tsec[kk]->GetTotalEnergy() ;
687  }
688  }
689 */
690  // Loop over hits retrieved from the library
691  for (unsigned int i=0; i<hits.getNhit(); i++) {
692 
693  // Get nPhotoElectrons and set edepositEM / edepositHAD accordingly
694  double nPhotoElectrons = hits.getNphotons(i);
695  // Apply scaling
696  nPhotoElectrons *= scale ;
697  if(isEM) {
698  // edepositEM = nPhotoElectrons*GeV;
699  edepositEM = nPhotoElectrons;
700  edepositHAD = 0.;
701  } else if(isHAD) {
702  edepositEM = 0.;
703  edepositHAD = nPhotoElectrons;
704  // edepositHAD = nPhotoElectrons*GeV;
705  }
706 
707  // Get hit position and time
708  double time = hits.getTime(i);
709  // math::XYZPoint position = hits.getHitPosition(i);
710 
711  // Get hit detID
712  unsigned int unitID = hits.getDetID(i);
713 
714  // Make the detID "rotation" from one sector to another taking into account the
715  // sectors of the impinging particle (theTrack) and of the particle that produced
716  // the 'hits' retrieved from shower library
717  unsigned int rotatedUnitID = rotateUnitID(unitID , theTrack , hits);
718  currentID.setID(rotatedUnitID, time, primaryID, 0);
719  // currentID.setID(unitID, time, primaryID, 0);
720 
721  // check if it is in the same unit and timeslice as the previous one
722  if (currentID == previousID) {
724  } else {
725  if (!checkHit()) currentHit = createNewHit();
726  }
727  } // End of loop over hits
728 
729  //Now kill the current track
730  if (ok) {
731  theTrack->SetTrackStatus(fStopAndKill);
732 #ifdef debugLog
733  LogDebug("ForwardSim") << "CastorSD::getFromLibrary:"
734  << "\n \"theTrack\" with TrackID() = "
735  << theTrack->GetTrackID()
736  << " and with energy "
737  << theTrack->GetTotalEnergy()
738  << " has been set to be killed" ;
739 #endif
740  G4TrackVector tv = *(aStep->GetSecondary());
741  for (unsigned int kk=0; kk<tv.size(); kk++) {
742  if (tv[kk]->GetVolume() == preStepPoint->GetPhysicalVolume()) {
743  tv[kk]->SetTrackStatus(fStopAndKill);
744 #ifdef debugLog
745  LogDebug("ForwardSim") << "CastorSD::getFromLibrary:"
746  << "\n tv[" << kk << "]->GetTrackID() = "
747  << tv[kk]->GetTrackID()
748  << " with energy "
749  << tv[kk]->GetTotalEnergy()
750  << " has been set to be killed" ;
751 #endif
752  }
753  }
754  }
755 }
#define LogDebug(id)
float edepositEM
Definition: CaloSD.h:121
int i
Definition: DBlmapReader.cc:9
G4int emPDG
Definition: CaloSD.h:136
void updateHit(CaloG4Hit *)
Definition: CaloSD.cc:440
double non_compensation_factor
Definition: CastorSD.h:53
G4ThreeVector posGlobal
Definition: CaloSD.h:113
uint32_t rotateUnitID(uint32_t, G4Track *, const CastorShowerEvent &)
Definition: CastorSD.cc:550
unsigned int getDetID(int i)
G4bool checkHit()
Definition: CaloSD.cc:331
void resetForNewPrimary(const G4ThreeVector &, double)
Definition: CaloSD.cc:454
float edepositHAD
Definition: CaloSD.h:121
G4int epPDG
Definition: CaloSD.h:136
float getTime(int i)
G4int gammaPDG
Definition: CaloSD.h:136
int setTrackID(G4Step *)
Definition: CastorSD.cc:527
float getNphotons(int i)
CaloHitID previousID
Definition: CaloSD.h:117
CaloG4Hit * currentHit
Definition: CaloSD.h:128
G4Track * theTrack
Definition: CaloSD.h:118
void setID(uint32_t unitID, double timeSlice, int trackID, uint16_t depth=0)
Definition: CaloHitID.cc:44
G4StepPoint * preStepPoint
Definition: CaloSD.h:120
static const G4LogicalVolume * GetVolume(const std::string &name)
CaloHitID currentID
Definition: CaloSD.h:117
CastorShowerLibrary * showerLibrary
Definition: CastorSD.h:47
CastorShowerEvent getShowerHits(G4Step *, bool &)
CaloG4Hit * createNewHit()
Definition: CaloSD.cc:363
unsigned int getNhit()
void CastorSD::initRun ( )
protectedvirtual

Reimplemented from CaloSD.

Definition at line 84 of file CastorSD.cc.

References CastorShowerLibrary::initParticleTable(), showerLibrary, and useShowerLibrary.

84  {
85  if (useShowerLibrary) {
86  // showerLibrary = new CastorShowerLibrary(name, cpv, p);
87  G4ParticleTable * theParticleTable = G4ParticleTable::GetParticleTable();
88  showerLibrary->initParticleTable(theParticleTable);
89  edm::LogInfo("ForwardSim") << "CastorSD::initRun: Using Castor Shower Library \n";
90  }
91 }
void initParticleTable(G4ParticleTable *)
bool useShowerLibrary
Definition: CastorSD.h:51
CastorShowerLibrary * showerLibrary
Definition: CastorSD.h:47
uint32_t CastorSD::rotateUnitID ( uint32_t  unitID,
G4Track *  track,
const CastorShowerEvent shower 
)
private

Definition at line 550 of file CastorSD.cc.

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

Referenced by getFromLibrary().

550  {
551 // ==============================================================
552 //
553 // o Exploit Castor phi symmetry to return newUnitID for
554 // shower hits based on track phi coordinate
555 //
556 // ==============================================================
557 
558  // Get 'track' phi:
559  float trackPhi = track->GetPosition().phi();
560  if(trackPhi<0) trackPhi += 2*M_PI ;
561  // Get phi from primary that gave rise to SL 'shower':
562  float showerPhi = shower.getPrimPhi();
563  if(showerPhi<0) showerPhi += 2*M_PI ;
564  // Delta phi:
565 
566  // Find the OctSector for which 'track' and 'shower' belong
567  int trackOctSector = (int) ( trackPhi / (M_PI/4) ) ;
568  int showerOctSector = (int) ( showerPhi / (M_PI/4) ) ;
569 
570  uint32_t newUnitID;
571  uint32_t sec = ( ( unitID>>4 ) & 0xF ) ;
572  uint32_t complement = ( unitID & 0xFFFFFF0F ) ;
573 
574  // edm::LogInfo("ForwardSim") << "\n CastorSD::rotateUnitID: "
575  // << "\n unitID = " << unitID
576  // << "\n sec = " << sec
577  // << "\n complement = " << complement ;
578 
579  // Get 'track' z:
580  float trackZ = track->GetPosition().z();
581 
582  int aux ;
583  int dSec = 2*(trackOctSector - showerOctSector) ;
584  // if(trackZ<0) // Good for revision 1.8 of CastorNumberingScheme
585  if(trackZ>0) // Good for revision 1.9 of CastorNumberingScheme
586  {
587  int sec1 = sec-dSec;
588  // sec -= dSec ;
589  if(sec1<0) sec1 += 16;
590  if(sec1>15) sec1 -= 16;
591  sec = (uint32_t)(sec1);
592  } else {
593  if( dSec<0 ) sec += 16 ;
594  sec += dSec ;
595  aux = (int) (sec/16) ;
596  sec -= aux*16 ;
597  }
598  sec = sec<<4 ;
599  newUnitID = complement | sec ;
600 
601 #ifdef debugLog
602  if(newUnitID != unitID) {
603  LogDebug("ForwardSim") << "\n CastorSD::rotateUnitID: "
604  << "\n unitID = " << unitID
605  << "\n newUnitID = " << newUnitID ;
606  }
607 #endif
608 
609  return newUnitID ;
610 
611 }
#define LogDebug(id)
#define M_PI
Definition: BFit3D.cc:3
float getPrimPhi() const
uint32_t CastorSD::setDetUnitId ( G4Step *  step)
virtual

Implements CaloSD.

Definition at line 509 of file CastorSD.cc.

References CastorNumberingScheme::getUnitID(), and numberingScheme.

509  {
510  return (numberingScheme == 0 ? 0 : numberingScheme->getUnitID(aStep));
511 }
CastorNumberingScheme * numberingScheme
Definition: CastorSD.h:46
virtual uint32_t getUnitID(const G4Step *aStep) const
void CastorSD::setNumberingScheme ( CastorNumberingScheme scheme)

Definition at line 515 of file CastorSD.cc.

References numberingScheme.

Referenced by CastorSD().

515  {
516 
517  if (scheme != 0) {
518  edm::LogInfo("ForwardSim") << "CastorSD: updates numbering scheme for "
519  << GetName();
520  if (numberingScheme) delete numberingScheme;
521  numberingScheme = scheme;
522  }
523 }
CastorNumberingScheme * numberingScheme
Definition: CastorSD.h:46
int CastorSD::setTrackID ( G4Step *  aStep)
private

Definition at line 527 of file CastorSD.cc.

References TrackInformation::getIDonCaloSurface(), CaloSD::preStepPoint, CaloSD::previousID, CaloSD::resetForNewPrimary(), CaloSD::theTrack, and CaloHitID::trackID().

Referenced by getFromLibrary().

527  {
528 
529  theTrack = aStep->GetTrack();
530 
531  double etrack = preStepPoint->GetKineticEnergy();
532  TrackInformation * trkInfo = (TrackInformation *)(theTrack->GetUserInformation());
533  int primaryID = trkInfo->getIDonCaloSurface();
534  if (primaryID == 0) {
535 #ifdef debugLog
536  edm::LogWarning("ForwardSim") << "CastorSD: Problem with primaryID **** set by force "
537  << "to TkID **** " << theTrack->GetTrackID();
538 #endif
539  primaryID = theTrack->GetTrackID();
540  }
541 
542  if (primaryID != previousID.trackID())
543  resetForNewPrimary(preStepPoint->GetPosition(), etrack);
544 
545  return primaryID;
546 }
int getIDonCaloSurface() const
void resetForNewPrimary(const G4ThreeVector &, double)
Definition: CaloSD.cc:454
int trackID() const
Definition: CaloHitID.h:25
CaloHitID previousID
Definition: CaloSD.h:117
G4Track * theTrack
Definition: CaloSD.h:118
G4StepPoint * preStepPoint
Definition: CaloSD.h:120

Member Data Documentation

double CastorSD::energyThresholdSL
private

Definition at line 52 of file CastorSD.h.

Referenced by CastorSD(), and getEnergyDeposit().

G4LogicalVolume* CastorSD::lvC3EF
private

Definition at line 48 of file CastorSD.h.

Referenced by CastorSD(), and getEnergyDeposit().

G4LogicalVolume * CastorSD::lvC3HF
private

Definition at line 48 of file CastorSD.h.

Referenced by CastorSD(), and getEnergyDeposit().

G4LogicalVolume * CastorSD::lvC4EF
private

Definition at line 48 of file CastorSD.h.

Referenced by CastorSD(), and getEnergyDeposit().

G4LogicalVolume * CastorSD::lvC4HF
private

Definition at line 48 of file CastorSD.h.

Referenced by CastorSD(), and getEnergyDeposit().

G4LogicalVolume* CastorSD::lvCAST
private

Definition at line 49 of file CastorSD.h.

Referenced by CastorSD(), and getEnergyDeposit().

double CastorSD::non_compensation_factor
private

Definition at line 53 of file CastorSD.h.

Referenced by CastorSD(), and getFromLibrary().

CastorNumberingScheme* CastorSD::numberingScheme
private

Definition at line 46 of file CastorSD.h.

Referenced by setDetUnitId(), and setNumberingScheme().

CastorShowerLibrary* CastorSD::showerLibrary
private

Definition at line 47 of file CastorSD.h.

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

bool CastorSD::useShowerLibrary
private

Definition at line 51 of file CastorSD.h.

Referenced by CastorSD(), getEnergyDeposit(), initRun(), and ~CastorSD().