CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HCalSD.cc
Go to the documentation of this file.
1 // File: HCalSD.cc
3 // Description: Sensitive Detector class for calorimeters
5 
17 
20 
21 #include "G4LogicalVolumeStore.hh"
22 #include "G4LogicalVolume.hh"
23 #include "G4Step.hh"
24 #include "G4Track.hh"
25 #include "G4ParticleTable.hh"
26 #include "G4VProcess.hh"
27 
28 #include <iostream>
29 #include <fstream>
30 #include <iomanip>
31 
32 #include "G4SystemOfUnits.hh"
33 
34 //#define DebugLog
35 //#define plotDebug
36 
37 HCalSD::HCalSD(G4String name, const DDCompactView & cpv,
39  edm::ParameterSet const & p, const SimTrackManager* manager) :
40  CaloSD(name, cpv, clg, p, manager,
41  p.getParameter<edm::ParameterSet>("HCalSD").getParameter<int>("TimeSliceUnit"),
42  p.getParameter<edm::ParameterSet>("HCalSD").getParameter<bool>("IgnoreTrackID")),
43  numberingFromDDD(0), numberingScheme(0), showerLibrary(0), hfshower(0),
44  showerParam(0), showerPMT(0), showerBundle(0), m_HEDarkening(0),
45  m_HFDarkening(0) {
46 
47  //static SimpleConfigurable<bool> on1(false, "HCalSD:UseBirkLaw");
48  //static SimpleConfigurable<double> bk1(0.013, "HCalSD:BirkC1");
49  //static SimpleConfigurable<double> bk2(0.0568,"HCalSD:BirkC2");
50  //static SimpleConfigurable<double> bk3(1.75, "HCalSD:BirkC3");
51  // Values from NIM 80 (1970) 239-244: as implemented in Geant3
52  //static SimpleConfigurable<bool> on2(true,"HCalSD:UseShowerLibrary");
54  useBirk = m_HC.getParameter<bool>("UseBirkLaw");
55  birk1 = m_HC.getParameter<double>("BirkC1")*(g/(MeV*cm2));
56  birk2 = m_HC.getParameter<double>("BirkC2");
57  birk3 = m_HC.getParameter<double>("BirkC3");
58  useShowerLibrary = m_HC.getParameter<bool>("UseShowerLibrary");
59  useParam = m_HC.getParameter<bool>("UseParametrize");
60  testNumber = m_HC.getParameter<bool>("TestNumberingScheme");
61  usePMTHit = m_HC.getParameter<bool>("UsePMTHits");
62  betaThr = m_HC.getParameter<double>("BetaThreshold");
63  eminHitHB = m_HC.getParameter<double>("EminHitHB")*MeV;
64  eminHitHE = m_HC.getParameter<double>("EminHitHE")*MeV;
65  eminHitHO = m_HC.getParameter<double>("EminHitHO")*MeV;
66  eminHitHF = m_HC.getParameter<double>("EminHitHF")*MeV;
67  useFibreBundle = m_HC.getParameter<bool>("UseFibreBundleHits");
68  deliveredLumi = m_HC.getParameter<double>("DelivLuminosity");
69  bool ageingFlagHE= m_HC.getParameter<bool>("HEDarkening");
70  bool ageingFlagHF= m_HC.getParameter<bool>("HFDarkening");
71  useHF = m_HC.getUntrackedParameter<bool>("UseHF",true);
72  bool forTBH2 = m_HC.getUntrackedParameter<bool>("ForTBH2",false);
73  useLayerWt = m_HC.getUntrackedParameter<bool>("UseLayerWt",false);
74  std::string file = m_HC.getUntrackedParameter<std::string>("WtFile","None");
75  edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
76  applyFidCut = m_HF.getParameter<bool>("ApplyFiducialCut");
77 
78 #ifdef DebugLog
79  LogDebug("HcalSim") << "***************************************************"
80  << "\n"
81  << "* *"
82  << "\n"
83  << "* Constructing a HCalSD with name " << name << "\n"
84  << "* *"
85  << "\n"
86  << "***************************************************";
87 #endif
88  edm::LogInfo("HcalSim") << "HCalSD:: Use of HF code is set to " << useHF
89  << "\nUse of shower parametrization set to "
90  << useParam << "\nUse of shower library is set to "
91  << useShowerLibrary << "\nUse PMT Hit is set to "
92  << usePMTHit << " with beta Threshold "<< betaThr
93  << "\nUSe of FibreBundle Hit set to "<<useFibreBundle
94  << "\n Use of Birks law is set to "
95  << useBirk << " with three constants kB = "
96  << birk1 << ", C1 = " << birk2 << ", C2 = " << birk3;
97  edm::LogInfo("HcalSim") << "HCalSD:: Suppression Flag " << suppressHeavy
98  << " protons below " << kmaxProton << " MeV,"
99  << " neutrons below " << kmaxNeutron << " MeV and"
100  << " ions below " << kmaxIon << " MeV\n"
101  << " Threshold for storing hits in HB: "
102  << eminHitHB << " HE: " << eminHitHE << " HO: "
103  << eminHitHO << " HF: " << eminHitHF << "\n"
104  << "Delivered luminosity for Darkening "
105  << deliveredLumi << " Flag (HE) " << ageingFlagHE
106  << " Flag (HF) " << ageingFlagHF << "\n"
107  << "Application of Fiducial Cut " << applyFidCut;
108 
109  numberingFromDDD = new HcalNumberingFromDDD(name, cpv);
110  HcalNumberingScheme* scheme;
111  if (testNumber || forTBH2)
112  scheme = dynamic_cast<HcalNumberingScheme*>(new HcalTestNumberingScheme(forTBH2));
113  else
114  scheme = new HcalNumberingScheme();
115  setNumberingScheme(scheme);
116 
117  const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
118  std::vector<G4LogicalVolume *>::const_iterator lvcite;
119  G4LogicalVolume* lv;
120  std::string attribute, value;
121  if (useHF) {
122  if (useParam) {
123  showerParam = new HFShowerParam(name, cpv, p);
124  } else {
125  if (useShowerLibrary) showerLibrary = new HFShowerLibrary(name, cpv, p);
126  hfshower = new HFShower(name, cpv, p, 0);
127  }
128 
129  // HF volume names
130  attribute = "Volume";
131  value = "HF";
132  DDSpecificsFilter filter0;
133  DDValue ddv0(attribute, value, 0);
134  filter0.setCriteria(ddv0, DDSpecificsFilter::equals);
135  DDFilteredView fv0(cpv);
136  fv0.addFilter(filter0);
137  hfNames = getNames(fv0);
138  fv0.firstChild();
139  DDsvalues_type sv0(fv0.mergedSpecifics());
140  std::vector<double> temp = getDDDArray("Levels",sv0);
141  edm::LogInfo("HcalSim") << "HCalSD: Names to be tested for " << attribute
142  << " = " << value << " has " << hfNames.size()
143  << " elements";
144  for (unsigned int i=0; i < hfNames.size(); ++i) {
145  G4String namv = hfNames[i];
146  lv = 0;
147  for(lvcite=lvs->begin(); lvcite!=lvs->end(); lvcite++)
148  if((*lvcite)->GetName()==namv) {
149  lv = (*lvcite);
150  break;
151  }
152  hfLV.push_back(lv);
153  int level = static_cast<int>(temp[i]);
154  hfLevels.push_back(level);
155  edm::LogInfo("HcalSim") << "HCalSD: HF[" << i << "] = " << hfNames[i]
156  << " LV " << hfLV[i] << " at level "
157  << hfLevels[i];
158  }
159 
160  // HF Fibre volume names
161  value = "HFFibre";
162  DDSpecificsFilter filter1;
163  DDValue ddv1(attribute,value,0);
164  filter1.setCriteria(ddv1, DDSpecificsFilter::equals);
165  DDFilteredView fv1(cpv);
166  fv1.addFilter(filter1);
167  fibreNames = getNames(fv1);
168  edm::LogInfo("HcalSim") << "HCalSD: Names to be tested for " << attribute
169  << " = " << value << ":";
170  for (unsigned int i=0; i<fibreNames.size(); ++i) {
171  G4String namv = fibreNames[i];
172  lv = 0;
173  for (lvcite = lvs->begin(); lvcite != lvs->end(); ++lvcite) {
174  if ((*lvcite)->GetName() == namv) {
175  lv = (*lvcite);
176  break;
177  }
178  }
179  fibreLV.push_back(lv);
180  edm::LogInfo("HcalSim") << "HCalSD: (" << i << ") " << fibreNames[i]
181  << " LV " << fibreLV[i];
182  }
183 
184  // HF PMT volume names
185  value = "HFPMT";
186  DDSpecificsFilter filter3;
187  DDValue ddv3(attribute,value,0);
189  DDFilteredView fv3(cpv);
190  fv3.addFilter(filter3);
191  std::vector<G4String> pmtNames = getNames(fv3);
192  edm::LogInfo("HcalSim") << "HCalSD: Names to be tested for " << attribute
193  << " = " << value << " have " << pmtNames.size()
194  << " entries";
195  for (unsigned int i=0; i<pmtNames.size(); ++i) {
196  G4String namv = pmtNames[i];
197  lv = 0;
198  for (lvcite = lvs->begin(); lvcite != lvs->end(); ++lvcite)
199  if ((*lvcite)->GetName() == namv) {
200  lv = (*lvcite);
201  break;
202  }
203  pmtLV.push_back(lv);
204  edm::LogInfo("HcalSim") << "HCalSD: (" << i << ") " << pmtNames[i]
205  << " LV " << pmtLV[i];
206  }
207  if (pmtNames.size() > 0) showerPMT = new HFShowerPMT (name, cpv, p);
208 
209  // HF Fibre bundles
210  value = "HFFibreBundleStraight";
211  DDSpecificsFilter filter4;
212  DDValue ddv4(attribute,value,0);
214  DDFilteredView fv4(cpv);
215  fv4.addFilter(filter4);
216  std::vector<G4String> fibreNames = getNames(fv4);
217  edm::LogInfo("HcalSim") << "HCalSD: Names to be tested for " << attribute
218  << " = " << value << " have " << fibreNames.size()
219  << " entries";
220  for (unsigned int i=0; i<fibreNames.size(); ++i) {
221  G4String namv = fibreNames[i];
222  lv = 0;
223  for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++)
224  if ((*lvcite)->GetName() == namv) {
225  lv = (*lvcite);
226  break;
227  }
228  fibre1LV.push_back(lv);
229  edm::LogInfo("HcalSim") << "HCalSD: (" << i << ") " << fibreNames[i]
230  << " LV " << fibre1LV[i];
231  }
232 
233  // Geometry parameters for HF
234  value = "HFFibreBundleConical";
235  DDSpecificsFilter filter5;
236  DDValue ddv5(attribute,value,0);
238  DDFilteredView fv5(cpv);
239  fv5.addFilter(filter5);
240  fibreNames = getNames(fv5);
241  edm::LogInfo("HcalSim") << "HCalSD: Names to be tested for " << attribute
242  << " = " << value << " have " << fibreNames.size()
243  << " entries";
244  for (unsigned int i=0; i<fibreNames.size(); ++i) {
245  G4String namv = fibreNames[i];
246  lv = 0;
247  for (lvcite = lvs->begin(); lvcite != lvs->end(); ++lvcite)
248  if ((*lvcite)->GetName() == namv) {
249  lv = (*lvcite);
250  break;
251  }
252  fibre2LV.push_back(lv);
253  edm::LogInfo("HcalSim") << "HCalSD: (" << i << ") " << fibreNames[i]
254  << " LV " << fibre2LV[i];
255  }
256  if (fibre1LV.size() > 0 || fibre2LV.size() > 0)
257  showerBundle = new HFShowerFibreBundle (name, cpv, p);
258 
259  attribute = "ReadOutName";
260  value = name;
261  DDSpecificsFilter filter6;
262  DDValue ddv6(attribute,value,0);
264  DDFilteredView fv6(cpv);
265  fv6.addFilter(filter6);
266  if (fv6.firstChild()) {
268  //Special Geometry parameters
269  gpar = getDDDArray("gparHF",sv);
270  edm::LogInfo("HFShower") << "HFShowerParam: " << gpar.size()
271  << " gpar (cm)";
272  for (unsigned int ig=0; ig<gpar.size(); ig++)
273  edm::LogInfo("HFShower") << "HFShowerParam: gpar[" << ig << "] = "
274  << gpar[ig]/cm << " cm";
275  } else {
276  edm::LogWarning("HFShower") << "HFShowerParam: cannot get filtered "
277  << " view for " << attribute << " matching "
278  << name;
279  }
280  }
281 
282  //Material list for HB/HE/HO sensitive detectors
283  attribute = "ReadOutName";
284  DDSpecificsFilter filter2;
285  DDValue ddv2(attribute,name,0);
287  DDFilteredView fv2(cpv);
288  fv2.addFilter(filter2);
289  bool dodet = fv2.firstChild();
290 
292  //Layer0 Weight
293  layer0wt = getDDDArray("Layer0Wt",sv);
294  edm::LogInfo("HcalSim") << "HCalSD: " << layer0wt.size() << " Layer0Wt";
295  for (unsigned int it=0; it<layer0wt.size(); ++it)
296  edm::LogInfo("HcalSim") << "HCalSD: [" << it << "] " << layer0wt[it];
297 
298  const G4MaterialTable * matTab = G4Material::GetMaterialTable();
299  std::vector<G4Material*>::const_iterator matite;
300  while (dodet) {
301  const DDLogicalPart & log = fv2.logicalPart();
302  G4String namx = log.name().name();
303  if (!isItHF(namx) && !isItFibre(namx)) {
304  bool notIn = true;
305  for (unsigned int i=0; i<matNames.size(); ++i) {
306  if (!strcmp(matNames[i].c_str(),log.material().name().name().c_str())){
307  notIn = false;
308  break;
309  }
310  }
311  if (notIn) {
312  namx = log.material().name().name();
313  matNames.push_back(namx);
314  G4Material* mat;
315  for (matite = matTab->begin(); matite != matTab->end(); ++matite) {
316  if ((*matite)->GetName() == namx) {
317  mat = (*matite);
318  break;
319  }
320  }
321  materials.push_back(mat);
322  }
323  }
324  dodet = fv2.next();
325  }
326 
327  edm::LogInfo("HcalSim") << "HCalSD: Material names for " << attribute
328  << " = " << name << ":";
329  for (unsigned int i=0; i<matNames.size(); ++i)
330  edm::LogInfo("HcalSim") << "HCalSD: (" << i << ") " << matNames[i]
331  << " pointer " << materials[i];
332 
333  mumPDG = mupPDG = 0;
334 
335  if (useLayerWt) readWeightFromFile(file);
336 
337  for (int i=0; i<9; ++i) hit_[i] = time_[i]= dist_[i] = 0;
338  hzvem = hzvhad = 0;
339 
340  if (ageingFlagHE) m_HEDarkening = new HEDarkening();
341  if (ageingFlagHF) m_HFDarkening = new HFDarkening();
342 #ifdef plotDebug
344 
345  if ( tfile.isAvailable() ) {
346  static const char * const labels[] = {"HB", "HE", "HO", "HF Absorber", "HF PMT",
347  "HF Absorber Long", "HF Absorber Short",
348  "HF PMT Long", "HF PMT Short"};
349  TFileDirectory hcDir = tfile->mkdir("ProfileFromHCalSD");
350  char name[20], title[60];
351  for (int i=0; i<9; ++i) {
352  sprintf (title, "Hit energy in %s", labels[i]);
353  sprintf (name, "HCalSDHit%d", i);
354  hit_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
355  sprintf (title, "Energy (MeV)");
356  hit_[i]->GetXaxis()->SetTitle(title);
357  hit_[i]->GetYaxis()->SetTitle("Hits");
358  sprintf (title, "Time of the hit in %s", labels[i]);
359  sprintf (name, "HCalSDTime%d", i);
360  time_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
361  sprintf (title, "Time (ns)");
362  time_[i]->GetXaxis()->SetTitle(title);
363  time_[i]->GetYaxis()->SetTitle("Hits");
364  sprintf (title, "Longitudinal profile in %s", labels[i]);
365  sprintf (name, "HCalSDDist%d", i);
366  dist_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
367  sprintf (title, "Distance (mm)");
368  dist_[i]->GetXaxis()->SetTitle(title);
369  dist_[i]->GetYaxis()->SetTitle("Hits");
370  }
371  if (useHF && (!useParam)) {
372  hzvem = hcDir.make<TH1F>("hzvem", "Longitudinal Profile (EM Part)",330,0.0,1650.0);
373  hzvem->GetXaxis()->SetTitle("Longitudinal Profile (EM Part)");
374  hzvhad = hcDir.make<TH1F>("hzvhad","Longitudinal Profile (Had Part)",330,0.0,1650.0);
375  hzvhad->GetXaxis()->SetTitle("Longitudinal Profile (Hadronic Part)");
376  }
377  }
378 #endif
379 }
380 
382 
384  if (numberingScheme) delete numberingScheme;
385  if (showerLibrary) delete showerLibrary;
386  if (hfshower) delete hfshower;
387  if (showerParam) delete showerParam;
388  if (showerPMT) delete showerPMT;
389  if (showerBundle) delete showerBundle;
390  if (m_HEDarkening) delete m_HEDarkening;
391  if (m_HFDarkening) delete m_HFDarkening;
392 }
393 
394 bool HCalSD::ProcessHits(G4Step * aStep, G4TouchableHistory * ) {
395 
396  NaNTrap( aStep ) ;
397 
398  if (aStep == NULL) {
399  return true;
400  } else {
401  G4LogicalVolume* lv =
402  aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
403  G4String nameVolume = lv->GetName();
404  if (isItHF(aStep)) {
405  G4int parCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
406  double weight(1.0);
407  if (m_HFDarkening != 0) {
408  G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
409  double r = hitPoint.perp()/CLHEP::cm;
410  double z = std::abs(hitPoint.z())/CLHEP::cm;
411  float dose_acquired = 0.;
412  if (z>=1100 && z <= 1300) {
413  int hfZLayer = (int)((z - 1100)/20);
414  if (hfZLayer > 9) hfZLayer = 9;
415  float normalized_lumi = m_HFDarkening->int_lumi(deliveredLumi);
416  for (int i = hfZLayer; i <= 9; ++i) {
417  dose_acquired = m_HFDarkening->dose(i,r);
418  weight *= m_HFDarkening->degradation(normalized_lumi*dose_acquired);
419  }
420  }
421 #ifdef DebugLog
422  LogDebug("HcalSim") << "HCalSD: HFLumiDarkening at r = " << r
423  << ", z = " << z << " Dose " << dose_acquired
424  << " weight " << weight;
425 #endif
426  }
427  if (useParam) {
428 #ifdef DebugLog
429  LogDebug("HcalSim") << "HCalSD: " << getNumberOfHits()
430  << " hits from parametrization in " << nameVolume
431  << " for Track " << aStep->GetTrack()->GetTrackID()
432  <<" (" << aStep->GetTrack()->GetDefinition()->GetParticleName()
433  <<")";
434 #endif
435  getFromParam(aStep, weight);
436 #ifdef DebugLog
437  LogDebug("HcalSim") << "HCalSD: " << getNumberOfHits()
438  << " hits afterParamS*";
439 #endif
440  } else {
441  bool notaMuon = true;
442  if (parCode == mupPDG || parCode == mumPDG ) notaMuon = false;
443  if (useShowerLibrary && notaMuon) {
444 #ifdef DebugLog
445  LogDebug("HcalSim") << "HCalSD: Starts shower library from "
446  << nameVolume << " for Track "
447  << aStep->GetTrack()->GetTrackID() << " ("
448  << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
449 #endif
450  getFromLibrary(aStep, weight);
451  } else if (isItFibre(lv)) {
452 #ifdef DebugLog
453  LogDebug("HcalSim") << "HCalSD: Hit at Fibre in " << nameVolume
454  << " for Track "
455  << aStep->GetTrack()->GetTrackID() <<" ("
456  << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
457 #endif
458  hitForFibre(aStep, weight);
459  }
460  }
461  } else if (isItPMT(lv)) {
462 #ifdef DebugLog
463  LogDebug("HcalSim") << "HCalSD: Hit from PMT parametrization from "
464  << nameVolume << " for Track "
465  << aStep->GetTrack()->GetTrackID() << " ("
466  << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
467 #endif
468  if (usePMTHit && showerPMT) getHitPMT(aStep);
469  } else if (isItStraightBundle(lv) || isItConicalBundle(lv)) {
470 #ifdef DebugLog
471  LogDebug("HcalSim") << "HCalSD: Hit from FibreBundle from "
472  << nameVolume << " for Track "
473  << aStep->GetTrack()->GetTrackID() << " ("
474  << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
475 #endif
478  } else {
479 #ifdef DebugLog
480  LogDebug("HcalSim") << "HCalSD: Hit from standard path from "
481  << nameVolume << " for Track "
482  << aStep->GetTrack()->GetTrackID() << " ("
483  << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
484 #endif
485  if (getStepInfo(aStep)) {
486 #ifdef plotDebug
487  if (edepositEM+edepositHAD > 0)
488  plotProfile(aStep, aStep->GetPreStepPoint()->GetPosition(),
489  edepositEM+edepositHAD,aStep->GetPostStepPoint()->GetGlobalTime(),0);
490 #endif
491  if (hitExists() == false && edepositEM+edepositHAD>0.) currentHit = createNewHit();
492  }
493  }
494  return true;
495  }
496 }
497 
498 double HCalSD::getEnergyDeposit(G4Step* aStep) {
499  double destep = aStep->GetTotalEnergyDeposit();
500  double weight = 1;
501  G4Track* theTrack = aStep->GetTrack();
502 
503  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
504  int depth = (touch->GetReplicaNumber(0))%10;
505  int det = ((touch->GetReplicaNumber(1))/1000)-3;
506  if (depth==0 && (det==0 || det==1)) weight = layer0wt[det];
507  if (useLayerWt) {
508  int lay = (touch->GetReplicaNumber(0)/10)%100 + 1;
509  G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
510  weight = layerWeight(det+3, hitPoint, depth, lay);
511  }
512 
513  if (m_HEDarkening !=0 && det == 1) {
514  int lay = (touch->GetReplicaNumber(0)/10)%100 + 1;
515  int ieta;
516  uint32_t detid = setDetUnitId(aStep);
517  if(testNumber) {
518  int det, z, depth, eta, phi, lay;
519  HcalTestNumbering::unpackHcalIndex(detid,det,z,depth,eta,phi,lay);
520  ieta = eta;
521  }
522  else ieta = HcalDetId(detid).ietaAbs();
523 #ifdef DebugLog
524  edm::LogInfo("HcalSimDark") << "HCalSD:HE_Darkening >>> ieta: "<< ieta //<< " vs. ietaAbs(): " << HcalDetId(detid).ietaAbs()
525  << " lay: " << lay-2;
526 #endif
527  float dweight = m_HEDarkening->degradation(deliveredLumi,ieta,lay-2);//NB:diff. layer count
528  weight *= dweight;
529 #ifdef DebugLog
530  edm::LogInfo("HcalSimDark") << "HCalSD: >>> Lumi: " << deliveredLumi
531  << " coefficient = " << dweight;
532 #endif
533  }
534 
535  if (suppressHeavy) {
536  TrackInformation * trkInfo = (TrackInformation *)(theTrack->GetUserInformation());
537  if (trkInfo) {
538  int pdg = theTrack->GetDefinition()->GetPDGEncoding();
539  if (!(trkInfo->isPrimary())) { // Only secondary particles
540  double ke = theTrack->GetKineticEnergy()/MeV;
541  if ( pdg/1000000000 == 1 && (pdg/10000)%100 > 0 &&
542  (pdg/10)%100 > 0 && ke <kmaxIon ) weight = 0;
543  if ((pdg == 2212) && (ke < kmaxProton)) weight = 0;
544  if ((pdg == 2112) && (ke < kmaxNeutron)) weight = 0;
545 #ifdef DebugLog
546  if (weight == 0)
547  edm::LogInfo("HcalSim") << "HCalSD:Ignore Track "
548  << theTrack->GetTrackID() << " Type "
549  << theTrack->GetDefinition()->GetParticleName()
550  << " Kinetic Energy " << ke << " MeV";
551 #endif
552  }
553  }
554  }
555 #ifdef DebugLog
556  double weight0 = weight;
557 #endif
558  if (useBirk) {
559  G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
560  if (isItScintillator(mat)) weight *= getAttenuation(aStep, birk1, birk2, birk3);
561  }
562  double wt1 = getResponseWt(theTrack);
563  double wt2 = theTrack->GetWeight();
564  /*
565  if(wt2 != 1.0) {
566  std::cout << "HCalSD: Detector " << det+3 << " Depth " << depth
567  << " weight= " << weight << " wt1= " << wt1
568  << " wt2= " << wt2 << std::endl;
569  std::cout << theTrack->GetDefinition()->GetParticleName()
570  << " " << theTrack->GetKineticEnergy()
571  << " Id=" << theTrack->GetTrackID()
572  << " IdP=" << theTrack->GetParentID();
573  const G4VProcess* pr = theTrack->GetCreatorProcess();
574  if(pr) std::cout << " from " << pr->GetProcessName();
575  std::cout << std::endl;
576  }
577  */
578 #ifdef DebugLog
579  edm::LogInfo("HcalSim") << "HCalSD: Detector " << det+3 << " Depth " << depth
580  << " weight " << weight0 << " " << weight << " " << wt1
581  << " " << wt2;
582 #endif
583  double edep = weight*wt1*destep;
584  if(wt2 > 0.0) { edep *= wt2; }
585  return edep;
586 }
587 
588 uint32_t HCalSD::setDetUnitId(G4Step * aStep) {
589 
590  G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
591  const G4VTouchable* touch = preStepPoint->GetTouchable();
592  G4ThreeVector hitPoint = preStepPoint->GetPosition();
593 
594  int depth = (touch->GetReplicaNumber(0))%10 + 1;
595  int lay = (touch->GetReplicaNumber(0)/10)%100 + 1;
596  int det = (touch->GetReplicaNumber(1))/1000;
597 
598  return setDetUnitId (det, hitPoint, depth, lay);
599 }
600 
602  if (scheme != 0) {
603  edm::LogInfo("HcalSim") << "HCalSD: updates numbering scheme for " << GetName();
604  if (numberingScheme) delete numberingScheme;
605  numberingScheme = scheme;
606  }
607 }
608 
610  G4ParticleTable * theParticleTable = G4ParticleTable::GetParticleTable();
611  G4String particleName;
612  mumPDG = theParticleTable->FindParticle(particleName="mu-")->GetPDGEncoding();
613  mupPDG = theParticleTable->FindParticle(particleName="mu+")->GetPDGEncoding();
614 #ifdef DebugLog
615  LogDebug("HcalSim") << "HCalSD: Particle code for mu- = " << mumPDG
616  << " for mu+ = " << mupPDG;
617 #endif
618  if (showerLibrary) showerLibrary->initRun(theParticleTable);
619  if (showerParam) showerParam->initRun(theParticleTable);
620  if (hfshower) hfshower->initRun(theParticleTable);
621 }
622 
623 bool HCalSD::filterHit(CaloG4Hit* aHit, double time) {
624  double threshold=0;
625  DetId theId(aHit->getUnitID());
626  switch (theId.subdetId()) {
627  case HcalBarrel:
628  threshold = eminHitHB; break;
629  case HcalEndcap:
630  threshold = eminHitHE; break;
631  case HcalOuter:
632  threshold = eminHitHO; break;
633  case HcalForward:
634  threshold = eminHitHF; break;
635  default:
636  break;
637  }
638  return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > threshold));
639 }
640 
641 
642 uint32_t HCalSD::setDetUnitId (int det, const G4ThreeVector& pos, int depth, int lay=1) {
643  uint32_t id = 0;
644  if (numberingFromDDD) {
645  //get the ID's as eta, phi, depth, ... indices
646  HcalNumberingFromDDD::HcalID tmp = numberingFromDDD->unitID(det, pos, depth, lay);
647  //get the ID
648  if (numberingScheme) id = numberingScheme->getUnitID(tmp);
649  }
650  return id;
651 }
652 
653 std::vector<double> HCalSD::getDDDArray(const std::string & str,
654  const DDsvalues_type & sv) {
655 #ifdef DebugLog
656  LogDebug("HcalSim") << "HCalSD:getDDDArray called for " << str;
657 #endif
658  DDValue value(str);
659  if (DDfetch(&sv,value)) {
660 #ifdef DebugLog
661  LogDebug("HcalSim") << value;
662 #endif
663  const std::vector<double> & fvec = value.doubles();
664  int nval = fvec.size();
665  if (nval < 1) {
666  edm::LogError("HcalSim") << "HCalSD : # of " << str << " bins " << nval
667  << " < 2 ==> illegal";
668  throw cms::Exception("Unknown", "HCalSD") << "nval < 2 for array " << str << "\n";
669  }
670 
671  return fvec;
672  } else {
673  edm::LogError("HcalSim") << "HCalSD : cannot get array " << str;
674  throw cms::Exception("Unknown", "HCalSD") << "cannot get array " << str << "\n";
675  }
676 }
677 
678 std::vector<G4String> HCalSD::getNames(DDFilteredView& fv) {
679 
680  std::vector<G4String> tmp;
681  bool dodet = fv.firstChild();
682  while (dodet) {
683  const DDLogicalPart & log = fv.logicalPart();
684  bool ok = true;
685 
686  for (unsigned int i=0; i<tmp.size(); ++i) {
687  if (!strcmp(tmp[i].c_str(), log.name().name().c_str())) {
688  ok = false;
689  break;
690  }
691  }
692  if (ok) tmp.push_back(log.name().name());
693  dodet = fv.next();
694  }
695  return tmp;
696 }
697 
698 bool HCalSD::isItHF(G4Step * aStep) {
699  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
700  int levels = (touch->GetHistoryDepth()) + 1;
701  for (unsigned int it=0; it < hfNames.size(); ++it) {
702  if (levels >= hfLevels[it]) {
703  G4LogicalVolume* lv = touch->GetVolume(levels-hfLevels[it])->GetLogicalVolume();
704  if (lv == hfLV[it]) return true;
705  }
706  }
707  return false;
708 }
709 
710 bool HCalSD::isItHF (G4String name) {
711  std::vector<G4String>::const_iterator it = hfNames.begin();
712  for (; it != hfNames.end(); ++it) if (name == *it) return true;
713  return false;
714 }
715 
716 bool HCalSD::isItFibre (G4LogicalVolume* lv) {
717  std::vector<G4LogicalVolume*>::const_iterator ite = fibreLV.begin();
718  for (; ite != fibreLV.end(); ++ite) if (lv == *ite) return true;
719  return false;
720 }
721 
722 bool HCalSD::isItFibre (G4String name) {
723  std::vector<G4String>::const_iterator it = fibreNames.begin();
724  for (; it != fibreNames.end(); ++it) if (name == *it) return true;
725  return false;
726 }
727 
728 bool HCalSD::isItPMT (G4LogicalVolume* lv) {
729  std::vector<G4LogicalVolume*>::const_iterator ite = pmtLV.begin();
730  for (; ite != pmtLV.end(); ++ite) if (lv == *ite) return true;
731  return false;
732 }
733 
734 bool HCalSD::isItStraightBundle (G4LogicalVolume* lv) {
735  std::vector<G4LogicalVolume*>::const_iterator ite = fibre1LV.begin();
736  for (; ite != fibre1LV.end(); ++ite) if (lv == *ite) return true;
737  return false;
738 }
739 
740 bool HCalSD::isItConicalBundle (G4LogicalVolume* lv) {
741  std::vector<G4LogicalVolume*>::const_iterator ite = fibre2LV.begin();
742  for (; ite != fibre2LV.end(); ++ite) if (lv == *ite) return true;
743  return false;
744 }
745 
746 bool HCalSD::isItScintillator (G4Material* mat) {
747  std::vector<G4Material*>::const_iterator ite = materials.begin();
748  for (; ite != materials.end(); ++ite) if (mat == *ite) return true;
749  return false;
750 }
751 
752 bool HCalSD::isItinFidVolume (G4ThreeVector& hitPoint) {
753  bool flag = true;
754  if (applyFidCut) {
755  int npmt = HFFibreFiducial:: PMTNumber(hitPoint);
756 #ifdef DebugLog
757  edm::LogInfo("HcalSim") << "HCalSD::isItinFidVolume:#PMT= " << npmt
758  << " for hit point " << hitPoint;
759 #endif
760  if (npmt <= 0) flag = false;
761  }
762 #ifdef DebugLog
763  edm::LogInfo("HcalSim") << "HCalSD::isItinFidVolume: point " << hitPoint
764  << " return flag " << flag;
765 #endif
766  return flag;
767 }
768 
769 void HCalSD::getFromLibrary (G4Step* aStep, double weight) {
770  preStepPoint = aStep->GetPreStepPoint();
771  theTrack = aStep->GetTrack();
772  int det = 5;
773  bool ok;
774 
775  std::vector<HFShowerLibrary::Hit> hits = showerLibrary->getHits(aStep, ok, weight, false);
776 
777  double etrack = preStepPoint->GetKineticEnergy();
778  int primaryID = setTrackID(aStep);
779 
780  // Reset entry point for new primary
781  posGlobal = preStepPoint->GetPosition();
782  resetForNewPrimary(posGlobal, etrack);
783 
784  G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
785  if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG) {
786  edepositEM = 1.*GeV;
787  edepositHAD = 0.;
788  } else {
789  edepositEM = 0.;
790  edepositHAD = 1.*GeV;
791  }
792 #ifdef DebugLog
793  edm::LogInfo("HcalSim") << "HCalSD::getFromLibrary " <<hits.size()
794  << " hits for " << GetName() << " of " << primaryID
795  << " with " << theTrack->GetDefinition()->GetParticleName()
796  << " of " << preStepPoint->GetKineticEnergy()/GeV << " GeV";
797 #endif
798  for (unsigned int i=0; i<hits.size(); ++i) {
799  G4ThreeVector hitPoint = hits[i].position;
800  if (isItinFidVolume (hitPoint)) {
801  int depth = hits[i].depth;
802  double time = hits[i].time;
803  unsigned int unitID = setDetUnitId(det, hitPoint, depth);
804  currentID.setID(unitID, time, primaryID, 0);
805 #ifdef plotDebug
806  plotProfile(aStep, hitPoint, 1.0*GeV, time, depth);
807  bool emType = false;
808  if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG)
809  emType = true;
810  plotHF(hitPoint,emType);
811 #endif
812 
813  // check if it is in the same unit and timeslice as the previous one
814  if (currentID == previousID) {
816  } else {
817  if (!checkHit()) currentHit = createNewHit();
818  }
819  }
820  }
821 
822  //Now kill the current track
823  if (ok) {
824  theTrack->SetTrackStatus(fStopAndKill);
825  G4TrackVector tv = *(aStep->GetSecondary());
826  for (unsigned int kk=0; kk<tv.size(); ++kk)
827  if (tv[kk]->GetVolume() == preStepPoint->GetPhysicalVolume())
828  tv[kk]->SetTrackStatus(fStopAndKill);
829  }
830 }
831 
832 void HCalSD::hitForFibre (G4Step* aStep, double weight) { // if not ParamShower
833 
834  preStepPoint = aStep->GetPreStepPoint();
835  theTrack = aStep->GetTrack();
836  int primaryID = setTrackID(aStep);
837 
838  int det = 5;
839  std::vector<HFShower::Hit> hits = hfshower->getHits(aStep, weight);
840 
841  G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
842  if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG) {
843  edepositEM = 1.*GeV;
844  edepositHAD = 0.;
845  } else {
846  edepositEM = 0.;
847  edepositHAD = 1.*GeV;
848  }
849 
850 #ifdef DebugLog
851  edm::LogInfo("HcalSim") << "HCalSD::hitForFibre " << hits.size()
852  << " hits for " << GetName() << " of " << primaryID
853  << " with " << theTrack->GetDefinition()->GetParticleName()
854  << " of " << preStepPoint->GetKineticEnergy()/GeV
855  << " GeV in detector type " << det;
856 #endif
857  if (hits.size() > 0) {
858  for (unsigned int i=0; i<hits.size(); ++i) {
859  G4ThreeVector hitPoint = hits[i].position;
860  if (isItinFidVolume (hitPoint)) {
861  int depth = hits[i].depth;
862  double time = hits[i].time;
863  unsigned int unitID = setDetUnitId(det, hitPoint, depth);
864  currentID.setID(unitID, time, primaryID, 0);
865 #ifdef plotDebug
866  plotProfile(aStep, hitPoint, edepositEM, time, depth);
867  bool emType = false;
868  if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG)
869  emType = true;
870  plotHF(hitPoint,emType);
871 #endif
872  // check if it is in the same unit and timeslice as the previous one
873  if (currentID == previousID) {
875  } else {
876  posGlobal = preStepPoint->GetPosition();
877  if (!checkHit()) currentHit = createNewHit();
878  }
879  }
880  }
881  }
882 }
883 
884 void HCalSD::getFromParam (G4Step* aStep, double weight) {
885  std::vector<HFShowerParam::Hit> hits = showerParam->getHits(aStep, weight);
886  int nHit = static_cast<int>(hits.size());
887 
888  if (nHit > 0) {
889  preStepPoint = aStep->GetPreStepPoint();
890  int primaryID = setTrackID(aStep);
891 
892  int det = 5;
893 #ifdef DebugLog
894  edm::LogInfo("HcalSim") << "HCalSD::getFromParam " << nHit << " hits for "
895  << GetName() << " of " << primaryID << " with "
896  << aStep->GetTrack()->GetDefinition()->GetParticleName()
897  << " of " << preStepPoint->GetKineticEnergy()/GeV
898  << " GeV in detector type " << det;
899 #endif
900  for (int i = 0; i<nHit; ++i) {
901  G4ThreeVector hitPoint = hits[i].position;
902  int depth = hits[i].depth;
903  double time = hits[i].time;
904  unsigned int unitID = setDetUnitId(det, hitPoint, depth);
905  currentID.setID(unitID, time, primaryID, 0);
906  edepositEM = hits[i].edep*GeV;
907  edepositHAD = 0.;
908 #ifdef plotDebug
909  plotProfile(aStep, hitPoint, edepositEM, time, depth);
910 #endif
911 
912  // check if it is in the same unit and timeslice as the previous one
913  if (currentID == previousID) {
915  } else {
916  posGlobal = preStepPoint->GetPosition();
917  if (!checkHit()) currentHit = createNewHit();
918  }
919  }
920  }
921 }
922 
923 void HCalSD::getHitPMT (G4Step * aStep) {
924 
925  preStepPoint = aStep->GetPreStepPoint();
926  theTrack = aStep->GetTrack();
927  double edep = showerPMT->getHits(aStep);
928 
929  if (edep >= 0) {
930  double etrack = preStepPoint->GetKineticEnergy();
931  int primaryID = 0;
932  if (etrack >= energyCut) {
933  primaryID = theTrack->GetTrackID();
934  } else {
935  primaryID = theTrack->GetParentID();
936  if (primaryID == 0) primaryID = theTrack->GetTrackID();
937  }
938  // Reset entry point for new primary
939  posGlobal = preStepPoint->GetPosition();
940  resetForNewPrimary(posGlobal, etrack);
941 
942  //
943  int det = static_cast<int>(HcalForward);
944  G4ThreeVector hitPoint = preStepPoint->GetPosition();
945  double rr = (hitPoint.x()*hitPoint.x() + hitPoint.y()*hitPoint.y());
946  double phi = (rr == 0. ? 0. :atan2(hitPoint.y(),hitPoint.x()));
947  double etaR = showerPMT->getRadius();
948  int depth = 1;
949  if (etaR < 0) {
950  depth = 2;
951  etaR =-etaR;
952  }
953  if (hitPoint.z() < 0) etaR =-etaR;
954 #ifdef DebugLog
955  edm::LogInfo("HcalSim") << "HCalSD::Hit for Detector " << det << " etaR "
956  << etaR << " phi " << phi/deg << " depth " <<depth;
957 #endif
958  double time = (aStep->GetPostStepPoint()->GetGlobalTime());
959  uint32_t unitID = 0;
960  if (numberingFromDDD) {
962  depth,1);
963  if (numberingScheme) unitID = numberingScheme->getUnitID(tmp);
964  }
965  currentID.setID(unitID, time, primaryID, 1);
966 
967  edepositHAD = aStep->GetTotalEnergyDeposit();
968  edepositEM =-edepositHAD + (edep*GeV);
969 #ifdef plotDebug
970  plotProfile(aStep, hitPoint, edep*GeV, time, depth);
971 #endif
972 #ifdef DebugLog
973  double beta = preStepPoint->GetBeta();
974  LogDebug("HcalSim") << "HCalSD::getHitPMT 1 hit for " << GetName()
975  << " of " << primaryID << " with "
976  << theTrack->GetDefinition()->GetParticleName()
977  << " of " << preStepPoint->GetKineticEnergy()/GeV
978  << " GeV with velocity " << beta << " UnitID "
979  << std::hex << unitID << std::dec;
980 #endif
981  // check if it is in the same unit and timeslice as the previous one
982  if (currentID == previousID) {
984  } else {
985  if (!checkHit()) currentHit = createNewHit();
986  }
987  }
988 }
989 
990 void HCalSD::getHitFibreBundle (G4Step* aStep, bool type) {
991  preStepPoint = aStep->GetPreStepPoint();
992  theTrack = aStep->GetTrack();
993  double edep = showerBundle->getHits(aStep, type);
994 
995  if (edep >= 0) {
996  double etrack = preStepPoint->GetKineticEnergy();
997  int primaryID = 0;
998  if (etrack >= energyCut) {
999  primaryID = theTrack->GetTrackID();
1000  } else {
1001  primaryID = theTrack->GetParentID();
1002  if (primaryID == 0) primaryID = theTrack->GetTrackID();
1003  }
1004  // Reset entry point for new primary
1005  posGlobal = preStepPoint->GetPosition();
1006  resetForNewPrimary(posGlobal, etrack);
1007 
1008  //
1009  int det = static_cast<int>(HcalForward);
1010  G4ThreeVector hitPoint = preStepPoint->GetPosition();
1011  double rr = (hitPoint.x()*hitPoint.x() + hitPoint.y()*hitPoint.y());
1012  double phi = (rr == 0. ? 0. :atan2(hitPoint.y(),hitPoint.x()));
1013  double etaR = showerBundle->getRadius();
1014  int depth = 1;
1015  if (etaR < 0) {
1016  depth = 2;
1017  etaR =-etaR;
1018  }
1019  if (hitPoint.z() < 0) etaR =-etaR;
1020 #ifdef DebugLog
1021  LogDebug("HcalSim") << "HCalSD::Hit for Detector " << det << " etaR "
1022  << etaR << " phi " << phi/deg << " depth " <<depth;
1023 #endif
1024  double time = (aStep->GetPostStepPoint()->GetGlobalTime());
1025  uint32_t unitID = 0;
1026  if (numberingFromDDD) {
1027  HcalNumberingFromDDD::HcalID tmp = numberingFromDDD->unitID(det,etaR,phi,depth,1);
1028  if (numberingScheme) unitID = numberingScheme->getUnitID(tmp);
1029  }
1030  if (type) currentID.setID(unitID, time, primaryID, 3);
1031  else currentID.setID(unitID, time, primaryID, 2);
1032 
1033  edepositHAD = aStep->GetTotalEnergyDeposit();
1034  edepositEM =-edepositHAD + (edep*GeV);
1035 #ifdef plotDebug
1036  plotProfile(aStep, hitPoint, edep*GeV, time, depth);
1037 #endif
1038 #ifdef DebugLog
1039  double beta = preStepPoint->GetBeta();
1040  LogDebug("HcalSim") << "HCalSD::getHitFibreBundle 1 hit for " << GetName()
1041  << " of " << primaryID << " with "
1042  << theTrack->GetDefinition()->GetParticleName()
1043  << " of " << preStepPoint->GetKineticEnergy()/GeV
1044  << " GeV with velocity " << beta << " UnitID "
1045  << std::hex << unitID << std::dec;
1046 #endif
1047  // check if it is in the same unit and timeslice as the previous one
1049  else if (!checkHit()) currentHit = createNewHit();
1050  } // non-zero energy deposit
1051 }
1052 
1053 int HCalSD::setTrackID (G4Step* aStep) {
1054  theTrack = aStep->GetTrack();
1055 
1056  double etrack = preStepPoint->GetKineticEnergy();
1057  TrackInformation * trkInfo = (TrackInformation *)(theTrack->GetUserInformation());
1058  int primaryID = trkInfo->getIDonCaloSurface();
1059  if (primaryID == 0) {
1060 #ifdef DebugLog
1061  edm::LogInfo("HcalSim") << "HCalSD: Problem with primaryID **** set by "
1062  << "force to TkID **** " <<theTrack->GetTrackID();
1063 #endif
1064  primaryID = theTrack->GetTrackID();
1065  }
1066 
1067  if (primaryID != previousID.trackID())
1068  resetForNewPrimary(preStepPoint->GetPosition(), etrack);
1069 
1070  return primaryID;
1071 }
1072 
1074 
1075  std::ifstream infile;
1076  int entry=0;
1077  infile.open(fName.c_str(), std::ios::in);
1078  if (infile) {
1079  int det, zside, etaR, phi, lay;
1080  double wt;
1081  while (infile >> det >> zside >> etaR >> phi >> lay >> wt) {
1082  uint32_t id = HcalTestNumbering::packHcalIndex(det,zside,1,etaR,phi,lay);
1083  layerWeights.insert(std::pair<uint32_t,double>(id,wt));
1084  ++entry;
1085 #ifdef DebugLog
1086  edm::LogInfo("HcalSim") << "HCalSD::readWeightFromFile:Entry " << entry
1087  << " ID " << std::hex << id << std::dec << " ("
1088  << det << "/" << zside << "/1/" << etaR << "/"
1089  << phi << "/" << lay << ") Weight " << wt;
1090 #endif
1091  }
1092  infile.close();
1093  }
1094  edm::LogInfo("HcalSim") << "HCalSD::readWeightFromFile: reads " << entry
1095  << " weights from " << fName;
1096  if (entry <= 0) useLayerWt = false;
1097 }
1098 
1099 double HCalSD::layerWeight(int det, const G4ThreeVector& pos, int depth, int lay) {
1100 
1101  double wt = 1.;
1102  if (numberingFromDDD) {
1103  //get the ID's as eta, phi, depth, ... indices
1105  depth, lay);
1106  uint32_t id = HcalTestNumbering::packHcalIndex(tmp.subdet, tmp.zside, 1,
1107  tmp.etaR, tmp.phis,tmp.lay);
1108  std::map<uint32_t,double>::const_iterator ite = layerWeights.find(id);
1109  if (ite != layerWeights.end()) wt = ite->second;
1110 #ifdef DebugLog
1111  edm::LogInfo("HcalSim") << "HCalSD::layerWeight: ID " << std::hex << id
1112  << std::dec << " (" << tmp.subdet << "/"
1113  << tmp.zside << "/1/" << tmp.etaR << "/"
1114  << tmp.phis << "/" << tmp.lay << ") Weight " <<wt;
1115 #endif
1116  }
1117  return wt;
1118 }
1119 
1120 void HCalSD::plotProfile(G4Step* aStep,const G4ThreeVector& global, double edep,
1121  double time, int id) {
1122 
1123  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
1124  static G4String modName[8] = {"HEModule", "HVQF" , "HBModule", "MBAT",
1125  "MBBT" , "MBBTC", "MBBT_R1P", "MBBT_R1M"};
1126  G4ThreeVector local;
1127  bool found=false;
1128  double depth=-2000;
1129  int idx = 4;
1130  for (int n=0; n<touch->GetHistoryDepth(); ++n) {
1131  G4String name = touch->GetVolume(n)->GetName();
1132 #ifdef DebugLog
1133  LogDebug("HcalSim") << "plotProfile Depth " << n << " Name " << name;
1134 #endif
1135  for (unsigned int ii=0; ii<8; ++ii) {
1136  if (name == modName[ii]) {
1137  found = true;
1138  int dn = touch->GetHistoryDepth() - n;
1139  local = touch->GetHistory()->GetTransform(dn).TransformPoint(global);
1140  if (ii == 0) {depth = local.z() - 4006.5; idx = 1;}
1141  else if (ii == 1) {depth = local.z() + 825.0; idx = 3;}
1142  else if (ii == 2) {depth = local.x() - 1775.; idx = 0;}
1143  else {depth = local.y() + 15.; idx = 2;}
1144  break;
1145  }
1146  }
1147  if (found) break;
1148  }
1149  if (!found) depth = std::abs(global.z()) - 11500;
1150 #ifdef DebugLog
1151  LogDebug("HcalSim") << "plotProfile Found " << found << " Global " << global
1152  << " Local " << local << " depth " << depth << " ID "
1153  << id << " EDEP " << edep << " Time " << time;
1154 #endif
1155  if (hit_[idx] != 0) hit_[idx]->Fill(edep);
1156  if (time_[idx] != 0) time_[idx]->Fill(time,edep);
1157  if (dist_[idx] != 0) dist_[idx]->Fill(depth,edep);
1158  int jd = 2*idx + id - 7;
1159  if (jd >= 0 && jd < 4) {
1160  jd += 5;
1161  if (hit_[jd] != 0) hit_[jd]->Fill(edep);
1162  if (time_[jd] != 0) time_[jd]->Fill(time,edep);
1163  if (dist_[jd] != 0) dist_[jd]->Fill(depth,edep);
1164  }
1165 }
1166 
1167 void HCalSD::plotHF(G4ThreeVector& hitPoint, bool emType) {
1168  double zv = std::abs(hitPoint.z()) - gpar[4];
1169  if (emType) {
1170  if (hzvem != 0) hzvem->Fill(zv);
1171  } else {
1172  if (hzvhad != 0) hzvhad->Fill(zv);
1173  }
1174 }
#define LogDebug(id)
float edepositEM
Definition: CaloSD.h:121
const double beta
double energyCut
Definition: CaloSD.h:123
type
Definition: HCALResponse.h:21
T getParameter(std::string const &) const
std::vector< double > layer0wt
Definition: HCalSD.h:92
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
bool useParam
Definition: HCalSD.h:88
double eminHitHE
Definition: HCalSD.h:89
HFShowerParam * showerParam
Definition: HCalSD.h:81
TH1F * time_[9]
Definition: HCalSD.h:98
G4int emPDG
Definition: CaloSD.h:136
double getRadius()
Definition: HFShowerPMT.cc:150
float int_lumi(float intlumi)
Definition: HFDarkening.cc:147
const DDLogicalPart & logicalPart() const
The logical-part of the current node in the filtered-view.
const std::vector< double > & doubles() const
a reference to the double-valued values stored in the given instance of DDValue
Definition: DDValue.cc:134
double kmaxNeutron
Definition: CaloSD.h:134
void addFilter(const DDFilter &, log_op op=AND)
void updateHit(CaloG4Hit *)
Definition: CaloSD.cc:440
void hitForFibre(G4Step *step, double weight)
Definition: HCalSD.cc:832
const N & name() const
Definition: DDBase.h:82
void getFromParam(G4Step *step, double weight)
Definition: HCalSD.cc:884
int getIDonCaloSurface() const
bool useLayerWt
Definition: HCalSD.h:86
Definition: CaloSD.h:42
void initRun(G4ParticleTable *)
Definition: HFShower.cc:478
bool isItStraightBundle(G4LogicalVolume *)
Definition: HCalSD.cc:734
virtual bool ProcessHits(G4Step *, G4TouchableHistory *)
Definition: HCalSD.cc:394
std::vector< double > gpar
Definition: HCalSD.h:92
bool useFibreBundle
Definition: HCalSD.h:86
bool isItinFidVolume(G4ThreeVector &)
Definition: HCalSD.cc:752
double betaThr
Definition: HCalSD.h:87
double deliveredLumi
Definition: HCalSD.h:90
std::vector< Hit > getHits(G4Step *aStep, bool &ok, double weight, bool onlyLong=false)
virtual uint32_t getUnitID(const HcalNumberingFromDDD::HcalID &id)
HFShowerFibreBundle * showerBundle
Definition: HCalSD.h:83
std::vector< G4LogicalVolume * > fibre1LV
Definition: HCalSD.h:96
std::vector< G4LogicalVolume * > fibreLV
Definition: HCalSD.h:96
double eminHitHB
Definition: HCalSD.h:89
bool useShowerLibrary
Definition: HCalSD.h:88
std::vector< double > getDDDArray(const std::string &, const DDsvalues_type &)
Definition: HCalSD.cc:653
double getHits(G4Step *aStep, bool type)
#define NULL
Definition: scimark2.h:8
virtual uint32_t setDetUnitId(G4Step *step)
Definition: HCalSD.cc:588
double birk2
Definition: HCalSD.h:87
int setTrackID(G4Step *step)
Definition: HCalSD.cc:1053
void getFromLibrary(G4Step *step, double weight)
Definition: HCalSD.cc:769
bool usePMTHit
Definition: HCalSD.h:86
G4ThreeVector posGlobal
Definition: CaloSD.h:113
T eta() const
int ii
Definition: cuy.py:588
void setNumberingScheme(HcalNumberingScheme *)
Definition: HCalSD.cc:601
type of data representation of DDCompactView
Definition: DDCompactView.h:77
double birk1
Definition: HCalSD.h:87
double kmaxProton
Definition: CaloSD.h:134
bool DDfetch(const DDsvalues_type *, DDValue &)
helper for retrieving DDValues from DDsvalues_type *.
Definition: DDsvalues.cc:102
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
float float float z
static int PMTNumber(const G4ThreeVector &pe_effect)
bool isItHF(G4Step *)
Definition: HCalSD.cc:698
G4bool checkHit()
Definition: CaloSD.cc:331
static TrackerG4SimHitNumberingScheme & numberingScheme(const DDCompactView &cpv, const GeometricDet &det)
TH1F * hzvhad
Definition: HCalSD.h:98
double eminHitHF
Definition: HCalSD.h:89
bool isItConicalBundle(G4LogicalVolume *)
Definition: HCalSD.cc:740
static uint32_t packHcalIndex(int det, int z, int depth, int eta, int phi, int lay)
TH1F * hit_[9]
Definition: HCalSD.h:98
void resetForNewPrimary(const G4ThreeVector &, double)
Definition: CaloSD.cc:454
virtual void initRun()
Definition: HCalSD.cc:609
bool testNumber
Definition: HCalSD.h:86
std::vector< G4LogicalVolume * > pmtLV
Definition: HCalSD.h:96
std::vector< G4LogicalVolume * > hfLV
Definition: HCalSD.h:96
double kmaxIon
Definition: CaloSD.h:134
bool suppressHeavy
Definition: CaloSD.h:132
bool isItScintillator(G4Material *)
Definition: HCalSD.cc:746
bool next()
set current node to the next node in the filtered tree
HFShower * hfshower
Definition: HCalSD.h:80
bool useBirk
Definition: HCalSD.h:86
float edepositHAD
Definition: CaloSD.h:121
double birk3
Definition: HCalSD.h:87
int trackID() const
Definition: CaloHitID.h:25
virtual double getEnergyDeposit(G4Step *)
Definition: HCalSD.cc:498
G4int epPDG
Definition: CaloSD.h:136
void plotHF(G4ThreeVector &pos, bool emType)
Definition: HCalSD.cc:1167
std::vector< std::pair< unsigned int, DDValue > > DDsvalues_type
std::maps an index to a DDValue. The index corresponds to the index assigned to the name of the std::...
Definition: DDsvalues.h:19
float degradation(float mrad)
Definition: HFDarkening.cc:143
G4int gammaPDG
Definition: CaloSD.h:136
bool isAvailable() const
Definition: Service.h:46
bool applyFidCut
Definition: HCalSD.h:88
std::vector< Hit > getHits(G4Step *aStep, double weight)
Definition: HFShower.cc:71
TH1F * dist_[9]
Definition: HCalSD.h:98
virtual ~HCalSD()
Definition: HCalSD.cc:381
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
static void unpackHcalIndex(const uint32_t &idx, int &det, int &z, int &depth, int &eta, int &phi, int &lay)
A DDLogicalPart aggregates information concerning material, solid and sensitveness ...
Definition: DDLogicalPart.h:88
HcalNumberingScheme * numberingScheme
Definition: HCalSD.h:78
void getHitPMT(G4Step *step)
Definition: HCalSD.cc:923
T * make(const Args &...args) const
make new ROOT object
int getNumberOfHits()
Definition: CaloSD.cc:361
HFDarkening * m_HFDarkening
Definition: HCalSD.h:85
CaloHitID previousID
Definition: CaloSD.h:117
double getAttenuation(G4Step *aStep, double birk1, double birk2, double birk3)
Definition: CaloSD.cc:465
CaloG4Hit * currentHit
Definition: CaloSD.h:128
std::vector< int > hfLevels
Definition: HCalSD.h:93
virtual std::vector< std::string > getNames()
void plotProfile(G4Step *step, const G4ThreeVector &pos, double edep, double time, int id)
Definition: HCalSD.cc:1120
HEDarkening * m_HEDarkening
Definition: HCalSD.h:84
void NaNTrap(G4Step *step)
G4Track * theTrack
Definition: CaloSD.h:118
virtual G4bool getStepInfo(G4Step *aStep)
Definition: CaloSD.cc:250
bool isItPMT(G4LogicalVolume *)
Definition: HCalSD.cc:728
G4int mupPDG
Definition: HCalSD.h:91
double eminHitHO
Definition: HCalSD.h:89
double tmaxHit
Definition: CaloSD.h:123
int ietaAbs() const
get the absolute value of the cell ieta
Definition: HcalDetId.h:34
void setID(uint32_t unitID, double timeSlice, int trackID, uint16_t depth=0)
Definition: CaloHitID.cc:44
TH1F * hzvem
Definition: HCalSD.h:98
G4StepPoint * preStepPoint
Definition: CaloSD.h:120
Definition: DetId.h:18
static const G4LogicalVolume * GetVolume(const std::string &name)
TFileDirectory mkdir(const std::string &dir, const std::string &descr="")
create a new subdirectory
Definition: TFileService.h:69
CaloHitID currentID
Definition: CaloSD.h:117
std::vector< Hit > getHits(G4Step *aStep, double weight)
void readWeightFromFile(std::string)
Definition: HCalSD.cc:1073
tuple idx
DEBUGGING if hasattr(process,&quot;trackMonIterativeTracking2012&quot;): print &quot;trackMonIterativeTracking2012 D...
void getHitFibreBundle(G4Step *step, bool type)
Definition: HCalSD.cc:990
int ke
bool isPrimary() const
DDsvalues_type mergedSpecifics() const
double getHits(G4Step *aStep)
Definition: HFShowerPMT.cc:104
void initRun(G4ParticleTable *theParticleTable)
HcalNumberingFromDDD * numberingFromDDD
Definition: HCalSD.h:77
double layerWeight(int, const G4ThreeVector &, int, int)
Definition: HCalSD.cc:1099
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
float dose(int layer, float radius)
Definition: HFDarkening.cc:92
double getResponseWt(G4Track *)
Definition: CaloSD.cc:607
void initRun(G4ParticleTable *)
std::vector< G4String > fibreNames
Definition: HCalSD.h:94
list infile
Definition: EdgesToViz.py:90
G4bool hitExists()
Definition: CaloSD.cc:310
bool firstChild()
set the current node to the first child ...
std::vector< G4Material * > materials
Definition: HCalSD.h:95
virtual bool filterHit(CaloG4Hit *, double)
Definition: HCalSD.cc:623
void setCriteria(const DDValue &nameVal, comp_op, log_op l=AND, bool asString=true, bool merged=true)
Definition: DDFilter.cc:285
std::vector< G4String > hfNames
Definition: HCalSD.h:94
G4int mumPDG
Definition: HCalSD.h:91
std::map< uint32_t, double > layerWeights
Definition: HCalSD.h:97
tuple level
Definition: testEve_cfg.py:34
HFShowerPMT * showerPMT
Definition: HCalSD.h:82
uint32_t getUnitID() const
Definition: CaloG4Hit.h:69
int weight
Definition: histoStyle.py:50
HcalID unitID(int det, const CLHEP::Hep3Vector &pos, int depth, int lay=-1) const
float degradation(float intlumi, int ieta, int lay)
Definition: HEDarkening.cc:52
HCalSD(G4String, const DDCompactView &, SensitiveDetectorCatalog &, edm::ParameterSet const &, const SimTrackManager *)
Definition: HCalSD.cc:37
HFShowerLibrary * showerLibrary
Definition: HCalSD.h:79
std::vector< G4LogicalVolume * > fibre2LV
Definition: HCalSD.h:96
const std::string & name() const
Returns the name.
Definition: DDName.cc:87
const DDMaterial & material(void) const
Returns a reference object of the material this LogicalPart is made of.
CaloG4Hit * createNewHit()
Definition: CaloSD.cc:363
std::vector< G4String > matNames
Definition: HCalSD.h:94
bool isItFibre(G4LogicalVolume *)
Definition: HCalSD.cc:716
double getEnergyDeposit() const
Definition: CaloG4Hit.h:81
bool useHF
Definition: HCalSD.h:88
Definition: DDAxes.h:10
The DDGenericFilter is a runtime-parametrized Filter looking on DDSpecifcs.
Definition: DDFilter.h:37