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 (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 (G4ThreeVector, double)
 
G4ThreeVector setToLocal (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 *, 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 30 of file CastorSD.cc.

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

33  :
34  CaloSD(name, cpv, clg, p, manager), numberingScheme(0), lvC3EF(0),
35  lvC3HF(0), lvC4EF(0), lvC4HF(0), lvCAST(0) {
36 
37  edm::ParameterSet m_CastorSD = p.getParameter<edm::ParameterSet>("CastorSD");
38  useShowerLibrary = m_CastorSD.getParameter<bool>("useShowerLibrary");
39  energyThresholdSL = m_CastorSD.getParameter<double>("minEnergyInGeVforUsingSLibrary");
40  energyThresholdSL = energyThresholdSL*GeV; // Convert GeV => MeV
41 
42  non_compensation_factor = m_CastorSD.getParameter<double>("nonCompensationFactor");
43 
44  if (useShowerLibrary) showerLibrary = new CastorShowerLibrary(name, p);
45 
47 
48  edm::LogInfo("ForwardSim")
49  << "***************************************************\n"
50  << "* *\n"
51  << "* Constructing a CastorSD with name " << GetName() << "\n"
52  << "* *\n"
53  << "***************************************************";
54 
55  const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
56  std::vector<G4LogicalVolume*>::const_iterator lvcite;
57  for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
58  if (strcmp(((*lvcite)->GetName()).c_str(),"C3EF") == 0) lvC3EF = (*lvcite);
59  if (strcmp(((*lvcite)->GetName()).c_str(),"C3HF") == 0) lvC3HF = (*lvcite);
60  if (strcmp(((*lvcite)->GetName()).c_str(),"C4EF") == 0) lvC4EF = (*lvcite);
61  if (strcmp(((*lvcite)->GetName()).c_str(),"C4HF") == 0) lvC4HF = (*lvcite);
62  if (strcmp(((*lvcite)->GetName()).c_str(),"CAST") == 0) lvCAST = (*lvcite);
63  if (lvC3EF != 0 && lvC3HF != 0 && lvC4EF != 0 && lvC4HF != 0 && lvCAST != 0) break;
64  }
65  edm::LogInfo("ForwardSim") << "CastorSD:: LogicalVolume pointers\n"
66  << lvC3EF << " for C3EF; " << lvC3HF
67  << " for C3HF; " << lvC4EF << " for C4EF; "
68  << lvC4HF << " for C4HF; "
69  << lvCAST << " for CAST. " << std::endl;
70 
71  // if(useShowerLibrary) edm::LogInfo("ForwardSim") << "\n Using Castor Shower Library \n";
72 
73 }
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:21
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:514
G4LogicalVolume * lvC3EF
Definition: CastorSD.h:48
CastorSD::~CastorSD ( )
virtual

Definition at line 77 of file CastorSD.cc.

References showerLibrary, and useShowerLibrary.

77  {
78  if (useShowerLibrary) delete showerLibrary;
79 }
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 94 of file CastorSD.cc.

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

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

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

Reimplemented from CaloSD.

Definition at line 83 of file CastorSD.cc.

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

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

Definition at line 549 of file CastorSD.cc.

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

Referenced by getFromLibrary().

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

Implements CaloSD.

Definition at line 508 of file CastorSD.cc.

References CastorNumberingScheme::getUnitID(), and numberingScheme.

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

Definition at line 514 of file CastorSD.cc.

References numberingScheme.

Referenced by CastorSD().

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

Definition at line 526 of file CastorSD.cc.

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

Referenced by getFromLibrary().

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

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().