CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CastorSD.cc
Go to the documentation of this file.
1 // File: CastorSD.cc
3 // Date: 02.04
4 // UpDate: 07.04 - C3TF & C4TF semi-trapezoid added
5 // Description: Sensitive Detector class for Castor
7 
8 // Added by WC
9 #include "SimG4Core/Notification/interface/TrackInformation.h"
10 
11 #include "SimG4CMS/Forward/interface/CastorSD.h"
12 //#include "SimDataFormats/CaloHit/interface/CastorShowerEvent.h"
13 #include "FWCore/MessageLogger/interface/MessageLogger.h"
14 
15 #include "G4SDManager.hh"
16 #include "G4Step.hh"
17 #include "G4Track.hh"
18 #include "G4VProcess.hh"
19 
20 #include "G4ios.hh"
21 #include "G4Cerenkov.hh"
22 #include "G4LogicalVolumeStore.hh"
23 
24 #include "CLHEP/Units/GlobalSystemOfUnits.h"
25 #include "Randomize.hh"
26 #include "G4Poisson.hh"
27 
28 //#define debugLog
29 
30 CastorSD::CastorSD(G4String name, const DDCompactView & cpv,
31  SensitiveDetectorCatalog & clg,
32  edm::ParameterSet const & p,
33  const SimTrackManager* manager) :
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 
46  setNumberingScheme(new CastorNumberingScheme());
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 }
74 
75 //=============================================================================================
76 
77 CastorSD::~CastorSD() {
78  if (useShowerLibrary) delete showerLibrary;
79 }
80 
81 //=============================================================================================
82 
83 void CastorSD::initRun(){
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 }
91 
92 //=============================================================================================
93 
94 double CastorSD::getEnergyDeposit(G4Step * aStep) {
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  float effPMTandTransport = 0.19;
301  /* for test HF geometry volumes:
302  if(nameVolume == "GF2Q" || nameVolume == "GFNQ" ||
303  nameVolume == "GR2Q" || nameVolume == "GRNQ")
304  effPMTandTransport = 0.19;
305  */
306 
307  float thFullRefl = 23.; /* 23.dergee */
308  float thFullReflRad = thFullRefl*pi/180.;
309 
310  /* default for Castor nameVolume == "CASF" or (C3TF & C4TF) */
311  float thFibDir = 45.; /* .dergee */
312  /* for test HF geometry volumes:
313  if(nameVolume == "GF2Q" || nameVolume == "GFNQ" ||
314  nameVolume == "GR2Q" || nameVolume == "GRNQ")
315  thFibDir = 0.0; // .dergee
316  */
317  float thFibDirRad = thFibDir*pi/180.;
318  /* */
319  /* */
320 
321  // at which theta the point is located:
322  // float th1 = hitPoint.theta();
323 
324  // theta of charged particle in LabRF(hit momentum direction):
325  float costh =hit_mom.z()/sqrt(hit_mom.x()*hit_mom.x()+
326  hit_mom.y()*hit_mom.y()+
327  hit_mom.z()*hit_mom.z());
328  if (zint < 0) costh = -costh;
329  float th = acos(std::min(std::max(costh,float(-1.)),float(1.)));
330 
331  // just in case (can do bot use):
332  if (th < 0.) th += twopi;
333 
334 
335 
336  // theta of cone with Cherenkov photons w.r.t.direction of charged part.:
337  float costhcher =1./(nMedium*beta);
338  float thcher = acos(std::min(std::max(costhcher,float(-1.)),float(1.)));
339 
340  // diff thetas of charged part. and quartz direction in LabRF:
341  float DelFibPart = fabs(th - thFibDirRad);
342 
343  // define real distances:
344  float d = fabs(tan(th)-tan(thFibDirRad));
345 
346  // float a = fabs(tan(thFibDirRad)-tan(thFibDirRad+thFullReflRad));
347  // float r = fabs(tan(th)-tan(th+thcher));
348 
349  float a = tan(thFibDirRad)+tan(fabs(thFibDirRad-thFullReflRad));
350  float r = tan(th)+tan(fabs(th-thcher));
351 
352 
353  // define losses d_qz in cone of full reflection inside quartz direction
354  float d_qz;
355  float variant;
356 
357  if(DelFibPart > (thFullReflRad + thcher) ) {d_qz = 0.; variant=0.;}
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.; variant=1.;}
363  // if((DelFibPart + thcher) < thFullReflRad ) {d_qz = 1.; variant=1.;}
364  // if((d+r) < a ) {d_qz = 1.; variant=1.;}
365  else {
366  if((thFibDirRad + thFullReflRad) < (th + thcher) &&
367  (thFibDirRad - thFullReflRad) > (th - thcher) ) {
368  // if((thcher - DelFibPart ) > thFullReflRad )
369  // if((r-d ) > a )
370  d_qz = 0.; variant=2.;
371 
372  } else {
373  // if((thcher + DelFibPart ) > thFullReflRad &&
374  // thcher < (DelFibPart+thFullReflRad) )
375  // {
376  d_qz = 0.; variant=3.;
377 
378 
379  // use crossed length of circles(cone projection)
380  // dC1/dC2 :
381  float arg_arcos = 0.;
382  float tan_arcos = 2.*a*d;
383  if(tan_arcos != 0.) arg_arcos =(r*r-a*a-d*d)/tan_arcos;
384  arg_arcos = fabs(arg_arcos);
385  float th_arcos = acos(std::min(std::max(arg_arcos,float(-1.)),float(1.)));
386  d_qz = th_arcos/pi/2.;
387  d_qz = fabs(d_qz);
388 
389 
390 
391  // }
392  // else
393  // {
394  // d_qz = 0.; variant=4.;
395  //#ifdef debugLog
396  // std::cout <<" ===============>variant 4 information: <===== " <<std::endl;
397  // std::cout <<" !!!!!!!!!!!!!!!!!!!!!! variant = " << variant <<std::endl;
398  //#endif
399  //
400  // }
401  }
402  }
403  }
404 
405 
406  // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
407 
408  if(charge != 0. && beta > bThreshold && d_qz != 0. ) {
409 
410  meanNCherPhot = 370.*charge*charge*
411  ( 1. - 1./(nMedium*nMedium*beta*beta) )*
412  photEnSpectrDE*stepl;
413 
414  // dLamdX:
415  // meanNCherPhot = (2.*pi/137.)*charge*charge*
416  // ( 1. - 1./(nMedium*nMedium*beta*beta) )*
417  // photEnSpectrDL*stepl;
418 
419 
420  // NCherPhot = meanNCherPhot;
421  // Poisson:
422  G4int poissNCherPhot = (G4int) G4Poisson(meanNCherPhot);
423 
424  if(poissNCherPhot < 0) poissNCherPhot = 0;
425 
426  NCherPhot = poissNCherPhot * effPMTandTransport * d_qz;
427 
428 
429 
430 #ifdef debugLog
431  float thgrad = th*180./pi;
432  float thchergrad = thcher*180./pi;
433  float DelFibPartgrad = DelFibPart*180./pi;
434  LogDebug("ForwardSim") << " ==============================> start all "
435  << "information:<========= \n" << " =====> for "
436  << "test:<=== \n" << " variant = " << variant
437  << "\n thgrad = " << thgrad << "\n thchergrad "
438  << "= " << thchergrad << "\n DelFibPartgrad = "
439  << DelFibPartgrad << "\n d_qz = " << d_qz
440  << "\n =====> Start Step Information <=== \n"
441  << " ===> calo preStepPoint info <=== \n"
442  << " hitPoint = " << hitPoint << "\n"
443  << " hitMom = " << hit_mom << "\n"
444  << " stepControlFlag = " << stepControlFlag
445  // << "\n curprocess = " << curprocess << "\n"
446  // << " nameProcess = " << nameProcess
447  << "\n charge = " << charge << "\n"
448  << " beta = " << beta << "\n"
449  << " bThreshold = " << bThreshold << "\n"
450  << " thgrad =" << thgrad << "\n"
451  << " effPMTandTransport=" << effPMTandTransport
452  // << "\n volume = " << name
453  << "\n nameVolume = " << nameVolume << "\n"
454  << " nMedium = " << nMedium << "\n"
455  // << " rad length = " << rad << "\n"
456  // << " material = " << mat << "\n"
457  << " stepl = " << stepl << "\n"
458  << " photEnSpectrDE = " << photEnSpectrDE <<"\n"
459  << " edep = " << edep << "\n"
460  << " ===> calo theTrack info <=== " << "\n"
461  << " particleType = " << particleType << "\n"
462  << " primaryID = " << primaryID << "\n"
463  << " entot= " << theTrack->GetTotalEnergy() << "\n"
464  << " vert_eta= " << eta << "\n"
465  << " vert_phi= " << phi << "\n"
466  << " vert_mom= " << vert_mom << "\n"
467  << " ===> calo hit preStepPointinfo <=== "<<"\n"
468  << " local point = " << localPoint << "\n"
469  << " ==============================> final info"
470  << ": <=== \n"
471  << " meanNCherPhot = " << meanNCherPhot << "\n"
472  << " poissNCherPhot = " << poissNCherPhot <<"\n"
473  << " NCherPhot = " << NCherPhot;
474 #endif
475 
476  // Included by WC
477  // std::cout << "\n volume = " << name
478  // << "\n nameVolume = " << nameVolume << "\n"
479  // << "\n postvolume = " << postname
480  // << "\n postnameVolume = " << postnameVolume << "\n"
481  // << "\n particleType = " << particleType
482  // << "\n primaryID = " << primaryID << "\n";
483 
484  }
485  }
486 
487 
488 #ifdef debugLog
489  LogDebug("ForwardSim") << "CastorSD:: " << nameVolume
490  // << " Light Collection Efficiency " << weight
491  << " Weighted Energy Deposit " << edep/MeV
492  << " MeV\n";
493 #endif
494  // Temporary member for testing purpose only...
495  // unit_id = setDetUnitId(aStep);
496  // if(NCherPhot != 0) std::cout << "\n UnitID = " << unit_id << " ; NCherPhot = " << NCherPhot ;
497 
498  return NCherPhot;
499  }
500  }
501  return 0;
502 }
503 
504 //=======================================================================================
505 
506 uint32_t CastorSD::setDetUnitId(G4Step* aStep) {
507  return (numberingScheme == 0 ? 0 : numberingScheme->getUnitID(aStep));
508 }
509 
510 //=======================================================================================
511 
512 void CastorSD::setNumberingScheme(CastorNumberingScheme* scheme) {
513 
514  if (scheme != 0) {
515  edm::LogInfo("ForwardSim") << "CastorSD: updates numbering scheme for "
516  << GetName();
517  if (numberingScheme) delete numberingScheme;
518  numberingScheme = scheme;
519  }
520 }
521 
522 //=======================================================================================
523 
524 int CastorSD::setTrackID (G4Step* aStep) {
525 
526  theTrack = aStep->GetTrack();
527 
528  double etrack = preStepPoint->GetKineticEnergy();
529  TrackInformation * trkInfo = (TrackInformation *)(theTrack->GetUserInformation());
530  int primaryID = trkInfo->getIDonCaloSurface();
531  if (primaryID == 0) {
532 #ifdef debugLog
533  edm::LogWarning("ForwardSim") << "CastorSD: Problem with primaryID **** set by force "
534  << "to TkID **** " << theTrack->GetTrackID();
535 #endif
536  primaryID = theTrack->GetTrackID();
537  }
538 
539  if (primaryID != previousID.trackID())
540  resetForNewPrimary(preStepPoint->GetPosition(), etrack);
541 
542  return primaryID;
543 }
544 
545 //=======================================================================================
546 
547 uint32_t CastorSD::rotateUnitID(uint32_t unitID, G4Track* track, CastorShowerEvent shower) {
548 // ==============================================================
549 //
550 // o Exploit Castor phi symmetry to return newUnitID for
551 // shower hits based on track phi coordinate
552 //
553 // ==============================================================
554 
555  // Get 'track' phi:
556  float trackPhi = track->GetPosition().phi();
557  if(trackPhi<0) trackPhi += 2*M_PI ;
558  // Get phi from primary that gave rise to SL 'shower':
559  float showerPhi = shower.getPrimPhi();
560  if(showerPhi<0) showerPhi += 2*M_PI ;
561  // Delta phi:
562 
563  // Find the OctSector for which 'track' and 'shower' belong
564  int trackOctSector = (int) ( trackPhi / (M_PI/4) ) ;
565  int showerOctSector = (int) ( showerPhi / (M_PI/4) ) ;
566 
567  uint32_t newUnitID;
568  uint32_t sec = ( ( unitID>>4 ) & 0xF ) ;
569  uint32_t complement = ( unitID & 0xFFFFFF0F ) ;
570 
571  // edm::LogInfo("ForwardSim") << "\n CastorSD::rotateUnitID: "
572  // << "\n unitID = " << unitID
573  // << "\n sec = " << sec
574  // << "\n complement = " << complement ;
575 
576  // Get 'track' z:
577  float trackZ = track->GetPosition().z();
578 
579  int aux ;
580  int dSec = 2*(trackOctSector - showerOctSector) ;
581  // if(trackZ<0) // Good for revision 1.8 of CastorNumberingScheme
582  if(trackZ>0) // Good for revision 1.9 of CastorNumberingScheme
583  {
584  int sec1 = sec-dSec;
585  // sec -= dSec ;
586  if(sec1<0) sec1 += 16;
587  if(sec1>15) sec1 -= 16;
588  sec = (uint32_t)(sec1);
589  } else {
590  if( dSec<0 ) sec += 16 ;
591  sec += dSec ;
592  aux = (int) (sec/16) ;
593  sec -= aux*16 ;
594  }
595  sec = sec<<4 ;
596  newUnitID = complement | sec ;
597 
598 #ifdef debugLog
599  if(newUnitID != unitID) {
600  LogDebug("ForwardSim") << "\n CastorSD::rotateUnitID: "
601  << "\n unitID = " << unitID
602  << "\n newUnitID = " << newUnitID ;
603  }
604 #endif
605 
606  return newUnitID ;
607 
608 }
609 
610 //=======================================================================================
611 
612 void CastorSD::getFromLibrary (G4Step* aStep) {
613 
615 //
616 // Method to get hits from the Shower Library
617 //
618 // CastorShowerEvent hits returned by getShowerHits are used to
619 // replace the full simulation of the shower from theTrack
620 //
621 // "updateHit" save the Hits to a CaloG4Hit container
622 //
624 
625  preStepPoint = aStep->GetPreStepPoint();
626  theTrack = aStep->GetTrack();
627  bool ok;
628 
629  // **** Call method to retrieve hits from the ShowerLibrary ****
630  CastorShowerEvent hits = showerLibrary->getShowerHits(aStep, ok);
631 
632  double etrack = preStepPoint->GetKineticEnergy();
633  int primaryID = setTrackID(aStep);
634  // int primaryID = theTrack->GetTrackID();
635 
636  // Reset entry point for new primary
637  posGlobal = preStepPoint->GetPosition();
638  resetForNewPrimary(posGlobal, etrack);
639 
640  // Check whether track is EM or HAD
641  G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
642  bool isEM , isHAD ;
643  if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG) {
644  isEM = true ; isHAD = false;
645  } else {
646  isEM = false; isHAD = true ;
647  }
648 
649 #ifdef debugLog
650  // edm::LogInfo("ForwardSim") << "\n CastorSD::getFromLibrary: "
651  LogDebug("ForwardSim") << "\n CastorSD::getFromLibrary: "
652  << hits.getNhit() << " hits for " << GetName() << " from "
653  << theTrack->GetDefinition()->GetParticleName() << " of "
654  << preStepPoint->GetKineticEnergy()/GeV << " GeV and trackID "
655  << theTrack->GetTrackID() ;
656 #endif
657 
658  // Scale to correct energy
659  double E_track = preStepPoint->GetTotalEnergy() ;
660  double E_SLhit = hits.getPrimE() * GeV ;
661  double scale = E_track/E_SLhit ;
662 
663  //Non compensation
664  if (isHAD){
665  scale=scale*non_compensation_factor; // if hadronic extend the scale with the non-compensation factor
666  } else {
667  scale=scale; // if electromagnetic, don't do anything
668  }
669 
670 
671 /* double theTrackEnergy = theTrack->GetTotalEnergy() ;
672 
673  if(fabs(theTrackEnergy-E_track)>10.) {
674  edm::LogInfo("ForwardSim") << "\n TrackID = " << theTrack->GetTrackID()
675  << "\n theTrackEnergy = " << theTrackEnergy
676  << "\n preStepPointEnergy = " << E_track ;
677  G4TrackVector tsec = *(aStep->GetSecondary());
678  for (unsigned int kk=0; kk<tsec.size(); kk++) {
679  edm::LogInfo("ForwardSim") << "CastorSD::getFromLibrary:"
680  << "\n tsec[" << kk << "]->GetTrackID() = "
681  << tsec[kk]->GetTrackID()
682  << " with energy "
683  << tsec[kk]->GetTotalEnergy() ;
684  }
685  }
686 */
687  // Loop over hits retrieved from the library
688  for (unsigned int i=0; i<hits.getNhit(); i++) {
689 
690  // Get nPhotoElectrons and set edepositEM / edepositHAD accordingly
691  double nPhotoElectrons = hits.getNphotons(i);
692  // Apply scaling
693  nPhotoElectrons *= scale ;
694  if(isEM) {
695  // edepositEM = nPhotoElectrons*GeV;
696  edepositEM = nPhotoElectrons;
697  edepositHAD = 0.;
698  } else if(isHAD) {
699  edepositEM = 0.;
700  edepositHAD = nPhotoElectrons;
701  // edepositHAD = nPhotoElectrons*GeV;
702  }
703 
704  // Get hit position and time
705  double time = hits.getTime(i);
706  math::XYZPoint position = hits.getHitPosition(i);
707 
708  // Get hit detID
709  unsigned int unitID = hits.getDetID(i);
710 
711  // Make the detID "rotation" from one sector to another taking into account the
712  // sectors of the impinging particle (theTrack) and of the particle that produced
713  // the 'hits' retrieved from shower library
714  unsigned int rotatedUnitID = rotateUnitID(unitID , theTrack , hits);
715  currentID.setID(rotatedUnitID, time, primaryID, 0);
716  // currentID.setID(unitID, time, primaryID, 0);
717 
718  // check if it is in the same unit and timeslice as the previous one
719  if (currentID == previousID) {
720  updateHit(currentHit);
721  } else {
722  if (!checkHit()) currentHit = createNewHit();
723  }
724  } // End of loop over hits
725 
726  //Now kill the current track
727  if (ok) {
728  theTrack->SetTrackStatus(fStopAndKill);
729 #ifdef debugLog
730  LogDebug("ForwardSim") << "CastorSD::getFromLibrary:"
731  << "\n \"theTrack\" with TrackID() = "
732  << theTrack->GetTrackID()
733  << " and with energy "
734  << theTrack->GetTotalEnergy()
735  << " has been set to be killed" ;
736 #endif
737  G4TrackVector tv = *(aStep->GetSecondary());
738  for (unsigned int kk=0; kk<tv.size(); kk++) {
739  if (tv[kk]->GetVolume() == preStepPoint->GetPhysicalVolume()) {
740  tv[kk]->SetTrackStatus(fStopAndKill);
741 #ifdef debugLog
742  LogDebug("ForwardSim") << "CastorSD::getFromLibrary:"
743  << "\n tv[" << kk << "]->GetTrackID() = "
744  << tv[kk]->GetTrackID()
745  << " with energy "
746  << tv[kk]->GetTotalEnergy()
747  << " has been set to be killed" ;
748 #endif
749  }
750  }
751  }
752 }
753 
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:516
CaloSD::edepositEM
float edepositEM
Definition: CaloSD.h:119
beta
const double beta
Definition: CosmicMuonParameters.h:69
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
CastorShowerLibrary::initParticleTable
void initParticleTable(G4ParticleTable *)
Definition: CastorShowerLibrary.cc:162
i
int i
Definition: DBlmapReader.cc:9
CaloSD::emPDG
G4int emPDG
Definition: CaloSD.h:134
CaloSD::updateHit
void updateHit(CaloG4Hit *)
Definition: CaloSD.cc:453
CastorShowerEvent::getHitPosition
Point getHitPosition(int i)
Definition: CastorShowerEvent.h:55
TrackInformation::getIDonCaloSurface
int getIDonCaloSurface() const
Definition: TrackInformation.h:36
CaloSD
Definition: CaloSD.h:42
CastorSD::getEnergyDeposit
virtual double getEnergyDeposit(G4Step *)
Definition: CastorSD.cc:94
CastorSD::initRun
virtual void initRun()
Definition: CastorSD.cc:83
CastorSD::rotateUnitID
uint32_t rotateUnitID(uint32_t, G4Track *, CastorShowerEvent)
Definition: CastorSD.cc:547
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:161
abs
#define abs(x)
Definition: mlp_lapack.h:159
CastorSD::non_compensation_factor
double non_compensation_factor
Definition: CastorSD.h:53
NULL
#define NULL
Definition: scimark2.h:8
min
#define min(a, b)
Definition: mlp_lapack.h:161
CaloSD::posGlobal
G4ThreeVector posGlobal
Definition: CaloSD.h:111
CastorSD::numberingScheme
CastorNumberingScheme * numberingScheme
Definition: CastorSD.h:46
CastorSD::useShowerLibrary
bool useShowerLibrary
Definition: CastorSD.h:51
eta
T eta() const
Definition: Basic3DVectorLD.h:168
DDCompactView
type of data representation of DDCompactView
Definition: DDCompactView.h:81
CastorSD::CastorSD
CastorSD(G4String, const DDCompactView &, SensitiveDetectorCatalog &clg, edm::ParameterSet const &, const SimTrackManager *)
Definition: CastorSD.cc:30
DeDxDiscriminatorTools::charge
double charge(const std::vector< uint8_t > &Ampls)
Definition: DeDxDiscriminatorTools.cc:43
CastorSD::getFromLibrary
void getFromLibrary(G4Step *)
Definition: CastorSD.cc:612
CastorShowerEvent::getDetID
unsigned int getDetID(int i)
Definition: CastorShowerEvent.h:54
position
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
CaloSD::checkHit
G4bool checkHit()
Definition: CaloSD.cc:343
numberingScheme
static TrackerG4SimHitNumberingScheme & numberingScheme(const DDCompactView &cpv, const GeometricDet &det)
Definition: TkAccumulatingSensitiveDetector.cc:54
CastorSD::lvC4HF
G4LogicalVolume * lvC4HF
Definition: CastorSD.h:48
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
CaloSD::edepositHAD
float edepositHAD
Definition: CaloSD.h:119
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
CaloHitID::trackID
int trackID() const
Definition: CaloHitID.h:25
CaloSD::epPDG
G4int epPDG
Definition: CaloSD.h:134
CastorShowerEvent::getTime
float getTime(int i)
Definition: CastorShowerEvent.h:57
CaloSD::gammaPDG
G4int gammaPDG
Definition: CaloSD.h:134
CastorSD::setTrackID
int setTrackID(G4Step *)
Definition: CastorSD.cc:524
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
CastorShowerEvent::getNphotons
float getNphotons(int i)
Definition: CastorShowerEvent.h:56
CastorSD::lvC3HF
G4LogicalVolume * lvC3HF
Definition: CastorSD.h:48
CaloSD::previousID
CaloHitID previousID
Definition: CaloSD.h:115
CaloSD::currentHit
CaloG4Hit * currentHit
Definition: CaloSD.h:126
CaloSD::theTrack
G4Track * theTrack
Definition: CaloSD.h:116
CastorSD::~CastorSD
virtual ~CastorSD()
Definition: CastorSD.cc:77
CastorSD::energyThresholdSL
double energyThresholdSL
Definition: CastorSD.h:52
CaloHitID::setID
void setID(uint32_t unitID, double timeSlice, int trackID, uint16_t depth=0)
Definition: CaloHitID.cc:44
CaloSD::preStepPoint
G4StepPoint * preStepPoint
Definition: CaloSD.h:118
CastorSD::lvC4EF
G4LogicalVolume * lvC4EF
Definition: CastorSD.h:48
GetVolume
static const G4LogicalVolume * GetVolume(const std::string &name)
Definition: TrackingMaterialProducer.cc:33
M_PI
#define M_PI
Definition: BFit3D.cc:3
funct::log
Log< T >::type log(const T &t)
Definition: Log.h:22
CaloSD::currentID
CaloHitID currentID
Definition: CaloSD.h:115
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:13
CastorSD::lvCAST
G4LogicalVolume * lvCAST
Definition: CastorSD.h:49
CastorSD::showerLibrary
CastorShowerLibrary * showerLibrary
Definition: CastorSD.h:47
L1TEmulatorMonitor_cff.p
tuple p
Definition: L1TEmulatorMonitor_cff.py:48
CastorShowerLibrary::getShowerHits
CastorShowerEvent getShowerHits(G4Step *, bool &)
Definition: CastorShowerLibrary.cc:211
CastorSD::setDetUnitId
virtual uint32_t setDetUnitId(G4Step *step)
Definition: CastorSD.cc:506
CastorSD::setNumberingScheme
void setNumberingScheme(CastorNumberingScheme *scheme)
Definition: CastorSD.cc:512
CastorSD::lvC3EF
G4LogicalVolume * lvC3EF
Definition: CastorSD.h:48
dot
T dot(const Basic3DVector &v) const
Scalar product, or &quot;dot&quot; product, with a vector of same type.
Definition: Basic3DVectorLD.h:223
a
double a
Definition: hdecay.h:121
CastorNumberingScheme::getUnitID
virtual uint32_t getUnitID(const G4Step *aStep) const
Definition: CastorNumberingScheme.cc:54
CaloSD::resetForNewPrimary
void resetForNewPrimary(G4ThreeVector, double)
Definition: CaloSD.cc:468
pi
double pi
Definition: Vx3DHLTAnalyzer.h:63
CaloSD::createNewHit
CaloG4Hit * createNewHit()
Definition: CaloSD.cc:375
CastorShowerEvent::getNhit
unsigned int getNhit()
Definition: CastorShowerEvent.h:53
phi
Definition: DDAxes.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5