CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HFShowerLibrary.cc
Go to the documentation of this file.
1 // File: HFShowerLibrary.cc
3 // Description: Shower library for Very forward hadron calorimeter
5 
6 #include "SimG4CMS/Calo/interface/HFShowerLibrary.h"
7 #include "SimDataFormats/CaloHit/interface/HFShowerLibraryEventInfo.h"
8 #include "DetectorDescription/Core/interface/DDFilter.h"
9 #include "DetectorDescription/Core/interface/DDFilteredView.h"
10 #include "DetectorDescription/Core/interface/DDValue.h"
11 
12 #include "FWCore/Utilities/interface/Exception.h"
13 
14 #include "G4VPhysicalVolume.hh"
15 #include "G4NavigationHistory.hh"
16 #include "G4Step.hh"
17 #include "G4Track.hh"
18 #include "Randomize.hh"
19 #include "CLHEP/Units/GlobalSystemOfUnits.h"
20 
21 //#define DebugLog
22 
23 HFShowerLibrary::HFShowerLibrary(std::string & name, const DDCompactView & cpv,
24  edm::ParameterSet const & p) : fibre(0),hf(0),
25  emBranch(0),
26  hadBranch(0),
27  npe(0) {
28 
29 
30  edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
31  probMax = m_HF.getParameter<double>("ProbMax");
32 
33  edm::ParameterSet m_HS= p.getParameter<edm::ParameterSet>("HFShowerLibrary");
34  edm::FileInPath fp = m_HS.getParameter<edm::FileInPath>("FileName");
35  std::string pTreeName = fp.fullPath();
36  backProb = m_HS.getParameter<double>("BackProbability");
37  std::string emName = m_HS.getParameter<std::string>("TreeEMID");
38  std::string hadName = m_HS.getParameter<std::string>("TreeHadID");
39  std::string branchEvInfo = m_HS.getUntrackedParameter<std::string>("BranchEvt","HFShowerLibraryEventInfos_hfshowerlib_HFShowerLibraryEventInfo");
40  std::string branchPre = m_HS.getUntrackedParameter<std::string>("BranchPre","HFShowerPhotons_hfshowerlib_");
41  std::string branchPost = m_HS.getUntrackedParameter<std::string>("BranchPost","_R.obj");
42  verbose = m_HS.getUntrackedParameter<bool>("Verbosity",false);
43 
44  if (pTreeName.find(".") == 0) pTreeName.erase(0,2);
45  const char* nTree = pTreeName.c_str();
46  hf = TFile::Open(nTree);
47 
48  if (!hf->IsOpen()) {
49  edm::LogError("HFShower") << "HFShowerLibrary: opening " << nTree
50  << " failed";
51  throw cms::Exception("Unknown", "HFShowerLibrary")
52  << "Opening of " << pTreeName << " fails\n";
53  } else {
54  edm::LogInfo("HFShower") << "HFShowerLibrary: opening " << nTree
55  << " successfully";
56  }
57 
58  TTree* event = (TTree *) hf ->Get("Events");
59  if (event) {
60  std::string info = branchEvInfo + branchPost;
61  TBranch *evtInfo = event->GetBranch(info.c_str());
62  if (evtInfo) {
63  loadEventInfo(evtInfo);
64  } else {
65  edm::LogError("HFShower") << "HFShowerLibrary: HFShowerLibrayEventInfo"
66  << " Branch does not exist in Event";
67  throw cms::Exception("Unknown", "HFShowerLibrary")
68  << "Event information absent\n";
69  }
70  } else {
71  edm::LogError("HFShower") << "HFShowerLibrary: Events Tree does not "
72  << "exist";
73  throw cms::Exception("Unknown", "HFShowerLibrary")
74  << "Events tree absent\n";
75  }
76 
77  edm::LogInfo("HFShower") << "HFShowerLibrary: Library " << libVers
78  << " ListVersion " << listVersion
79  << " Events Total " << totEvents << " and "
80  << evtPerBin << " per bin";
81  edm::LogInfo("HFShower") << "HFShowerLibrary: Energies (GeV) with "
82  << nMomBin << " bins";
83  for (int i=0; i<nMomBin; i++)
84  edm::LogInfo("HFShower") << "HFShowerLibrary: pmom[" << i << "] = "
85  << pmom[i]/GeV << " GeV";
86 
87  std::string nameBr = branchPre + emName + branchPost;
88  emBranch = event->GetBranch(nameBr.c_str());
89  if (verbose) emBranch->Print();
90  nameBr = branchPre + hadName + branchPost;
91  hadBranch = event->GetBranch(nameBr.c_str());
92  if (verbose) hadBranch->Print();
93  edm::LogInfo("HFShower") << "HFShowerLibrary:Branch " << emName
94  << " has " << emBranch->GetEntries()
95  << " entries and Branch " << hadName
96  << " has " << hadBranch->GetEntries()
97  << " entries";
98  edm::LogInfo("HFShower") << "HFShowerLibrary::No packing information -"
99  << " Assume x, y, z are not in packed form";
100 
101  edm::LogInfo("HFShower") << "HFShowerLibrary: Maximum probability cut off "
102  << probMax << " Back propagation of light prob. "
103  << backProb ;
104 
105  G4String attribute = "ReadOutName";
106  G4String value = name;
107  DDSpecificsFilter filter;
108  DDValue ddv(attribute,value,0);
109  filter.setCriteria(ddv,DDSpecificsFilter::equals);
110  DDFilteredView fv(cpv);
111  fv.addFilter(filter);
112  bool dodet = fv.firstChild();
113  if (dodet) {
114  DDsvalues_type sv(fv.mergedSpecifics());
115 
116  //Radius (minimum and maximum)
117  int nR = -1;
118  std::vector<double> rTable = getDDDArray("rTable",sv,nR);
119  rMin = rTable[0];
120  rMax = rTable[nR-1];
121  edm::LogInfo("HFShower") << "HFShowerLibrary: rMIN " << rMin/cm
122  << " cm and rMax " << rMax/cm;
123 
124  //Delta phi
125  int nEta = -1;
126  std::vector<double> etaTable = getDDDArray("etaTable",sv,nEta);
127  int nPhi = nEta + nR - 2;
128  std::vector<double> phibin = getDDDArray("phibin",sv,nPhi);
129  dphi = phibin[nEta-1];
130  edm::LogInfo("HFShower") << "HFShowerLibrary: (Half) Phi Width of wedge "
131  << dphi/deg;
132 
133  //Special Geometry parameters
134  int ngpar = 7;
135  gpar = getDDDArray("gparHF",sv,ngpar);
136  edm::LogInfo("HFShower") << "HFShowerLibrary: " << ngpar << " gpar (cm)";
137  for (int ig=0; ig<ngpar; ig++)
138  edm::LogInfo("HFShower") << "HFShowerLibrary: gpar[" << ig << "] = "
139  << gpar[ig]/cm << " cm";
140  } else {
141  edm::LogError("HFShower") << "HFShowerLibrary: cannot get filtered "
142  << " view for " << attribute << " matching "
143  << name;
144  throw cms::Exception("Unknown", "HFShowerLibrary")
145  << "cannot match " << attribute << " to " << name <<"\n";
146  }
147 
148  fibre = new HFFibre(name, cpv, p);
149  emPDG = epPDG = gammaPDG = 0;
150  pi0PDG = etaPDG = nuePDG = numuPDG = nutauPDG= 0;
151  anuePDG= anumuPDG = anutauPDG = geantinoPDG = 0;
152 }
153 
154 HFShowerLibrary::~HFShowerLibrary() {
155  if (hf) hf->Close();
156  if (fibre) delete fibre; fibre = 0;
157 }
158 
159 void HFShowerLibrary::initRun(G4ParticleTable * theParticleTable) {
160 
161  G4String parName;
162  emPDG = theParticleTable->FindParticle(parName="e-")->GetPDGEncoding();
163  epPDG = theParticleTable->FindParticle(parName="e+")->GetPDGEncoding();
164  gammaPDG = theParticleTable->FindParticle(parName="gamma")->GetPDGEncoding();
165  pi0PDG = theParticleTable->FindParticle(parName="pi0")->GetPDGEncoding();
166  etaPDG = theParticleTable->FindParticle(parName="eta")->GetPDGEncoding();
167  nuePDG = theParticleTable->FindParticle(parName="nu_e")->GetPDGEncoding();
168  numuPDG = theParticleTable->FindParticle(parName="nu_mu")->GetPDGEncoding();
169  nutauPDG= theParticleTable->FindParticle(parName="nu_tau")->GetPDGEncoding();
170  anuePDG = theParticleTable->FindParticle(parName="anti_nu_e")->GetPDGEncoding();
171  anumuPDG= theParticleTable->FindParticle(parName="anti_nu_mu")->GetPDGEncoding();
172  anutauPDG= theParticleTable->FindParticle(parName="anti_nu_tau")->GetPDGEncoding();
173  geantinoPDG= theParticleTable->FindParticle(parName="geantino")->GetPDGEncoding();
174 #ifdef DebugLog
175  edm::LogInfo("HFShower") << "HFShowerLibrary: Particle codes for e- = "
176  << emPDG << ", e+ = " << epPDG << ", gamma = "
177  << gammaPDG << ", pi0 = " << pi0PDG << ", eta = "
178  << etaPDG << ", geantino = " << geantinoPDG
179  << "\n nu_e = " << nuePDG << ", nu_mu = "
180  << numuPDG << ", nu_tau = " << nutauPDG
181  << ", anti_nu_e = " << anuePDG << ", anti_nu_mu = "
182  << anumuPDG << ", anti_nu_tau = " << anutauPDG;
183 #endif
184 }
185 
186 
187 std::vector<HFShowerLibrary::Hit> HFShowerLibrary::getHits(G4Step * aStep,
188  bool & ok,
189  bool onlyLong) {
190 
191  G4StepPoint * preStepPoint = aStep->GetPreStepPoint();
192  G4StepPoint * postStepPoint = aStep->GetPostStepPoint();
193  G4Track * track = aStep->GetTrack();
194  // Get Z-direction
195  const G4DynamicParticle *aParticle = track->GetDynamicParticle();
196  G4ThreeVector momDir = aParticle->GetMomentumDirection();
197  // double mom = aParticle->GetTotalMomentum();
198 
199  G4ThreeVector hitPoint = preStepPoint->GetPosition();
200  G4String partType = track->GetDefinition()->GetParticleName();
201  int parCode = track->GetDefinition()->GetPDGEncoding();
202 
203  std::vector<HFShowerLibrary::Hit> hit;
204  ok = false;
205  if (parCode == pi0PDG || parCode == etaPDG || parCode == nuePDG ||
206  parCode == numuPDG || parCode == nutauPDG || parCode == anuePDG ||
207  parCode == anumuPDG || parCode == anutauPDG || parCode == geantinoPDG)
208  return hit;
209  ok = true;
210 
211  double tSlice = (postStepPoint->GetGlobalTime())/nanosecond;
212  double pin = preStepPoint->GetTotalEnergy();
213  double pz = momDir.z();
214  double zint = hitPoint.z();
215 
216  // if particle moves from interaction point or "backwards (halo)
217  int backward = 0;
218  if (pz * zint < 0.) backward = 1;
219 
220  double sphi = sin(momDir.phi());
221  double cphi = cos(momDir.phi());
222  double ctheta = cos(momDir.theta());
223  double stheta = sin(momDir.theta());
224 
225 #ifdef DebugLog
226  G4ThreeVector localPos = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
227  double zoff = localPos.z() + 0.5*gpar[1];
228  // if (zoff < 0) zoff = 0;
229  edm::LogInfo("HFShower") << "HFShowerLibrary: getHits " << partType
230  << " of energy " << pin/GeV << " GeV"
231  << " dir.orts " << momDir.x() << ", " <<momDir.y()
232  << ", " << momDir.z() << " Pos x,y,z = "
233  << hitPoint.x() << "," << hitPoint.y() << ","
234  << hitPoint.z() << " (" << zoff
235  << ") sphi,cphi,stheta,ctheta = " << sphi
236  << "," << cphi << ", " << stheta << "," << ctheta;
237 #endif
238 
239  if (parCode == emPDG || parCode == epPDG || parCode == gammaPDG ) {
240  if (pin<pmom[nMomBin-1]) {
241  interpolate(0, pin);
242  } else {
243  extrapolate(0, pin);
244  }
245  } else {
246  if (pin<pmom[nMomBin-1]) {
247  interpolate(1, pin);
248  } else {
249  extrapolate(1, pin);
250  }
251  }
252 
253  int nHit = 0;
254  HFShowerLibrary::Hit oneHit;
255  for (int i = 0; i < npe; i++) {
256  double zv = std::abs(pe[i].z()); // abs local z
257 #ifdef DebugLog
258  LogDebug("HFShower") << "HFShowerLibrary: Hit " << i << " " << pe[i] << " zv " << zv;
259 #endif
260  if (zv <= gpar[1] && pe[i].lambda() > 0 &&
261  (pe[i].z() >= 0 || (zv > gpar[0] && (!onlyLong)))) {
262  int depth = 1;
263  if (onlyLong) {
264  } else if (backward == 0) { // fully valid only for "front" particles
265  if (pe[i].z() < 0) depth = 2;// with "front"-simulated shower lib.
266  } else { // for "backward" particles - almost equal
267  double r = G4UniformRand(); // share between L and S fibers
268  if (r > 0.5) depth = 2;
269  }
270 
271 
272  // Updated coordinate transformation from local
273  // back to global using two Euler angles: phi and theta
274  double pex = pe[i].x();
275  double pey = pe[i].y();
276 
277  double xx = pex*ctheta*cphi - pey*sphi + zv*stheta*cphi;
278  double yy = pex*ctheta*sphi + pey*cphi + zv*stheta*sphi;
279  double zz = -pex*stheta + zv*ctheta;
280 
281  G4ThreeVector pos = hitPoint + G4ThreeVector(xx,yy,zz);
282  zv = std::abs(pos.z()) - gpar[4] - 0.5*gpar[1];
283  G4ThreeVector lpos = G4ThreeVector(pos.x(),pos.y(),zv);
284 
285  zv = fibre->zShift(lpos,depth,0); // distance to PMT !
286 
287  double r = pos.perp();
288  double p = fibre->attLength(pe[i].lambda());
289  double fi = pos.phi();
290  if (fi < 0) fi += twopi;
291  int isect = int(fi/dphi) + 1;
292  isect = (isect + 1) / 2;
293  double dfi = ((isect*2-1)*dphi - fi);
294  if (dfi < 0) dfi = -dfi;
295  double dfir = r * sin(dfi);
296 #ifdef DebugLog
297  LogDebug("HFShower") << "HFShowerLibrary: Position shift " << xx
298  << ", " << yy << ", " << zz << ": " << pos
299  << " R " << r << " Phi " << fi << " Section "
300  << isect << " R*Dfi " << dfir << " Dist " << zv;
301 #endif
302  zz = std::abs(pos.z());
303  double r1 = G4UniformRand();
304  double r2 = G4UniformRand();
305  double r3 = -9999.;
306  if(backward) r3 = G4UniformRand();
307 
308 #ifdef DebugLog
309  LogDebug("HFShower") << "HFShowerLibrary: rLimits " << rInside(r)
310  << " attenuation " << r1 <<":" << exp(-p*zv)
311  << " r2 " << r2 << " r3 " << r3 << " rDfi "
312  << gpar[5] << " zz "
313  << zz << " zLim " << gpar[4] << ":"
314  << gpar[4]+gpar[1] << "\n"
315  << " rInside(r) :" << rInside(r)
316  << " r1 <= exp(-p*zv) :" << (r1 <= exp(-p*zv))
317  << " r2 <= probMax :" << (r2 <= probMax)
318  << " r3 <= backProb :" << (r3 <= backProb)
319  << " dfir > gpar[5] :" << (dfir > gpar[5])
320  << " zz >= gpar[4] :" << (zz >= gpar[4])
321  << " zz <= gpar[4]+gpar[1] :"
322  << (zz <= gpar[4]+gpar[1]);
323 #endif
324  if (rInside(r) && r1 <= exp(-p*zv) && r2 <= probMax && dfir > gpar[5] &&
325  zz >= gpar[4] && zz <= gpar[4]+gpar[1] && r3 <= backProb &&
326  (depth != 2 || zz >= gpar[4]+gpar[0])) {
327  oneHit.position = pos;
328  oneHit.depth = depth;
329  oneHit.time = (tSlice+(pe[i].t())+(fibre->tShift(lpos,depth,1)));
330  hit.push_back(oneHit);
331 #ifdef DebugLog
332  LogDebug("HFShower") << "HFShowerLibrary: Final Hit " << nHit
333  <<" position " << (hit[nHit].position) <<" Depth "
334  << (hit[nHit].depth) <<" Time " << tSlice << ":"
335  << pe[i].t() << ":" <<fibre->tShift(lpos,depth,1)
336  << ":" << (hit[nHit].time);
337 #endif
338  nHit++;
339  }
340 #ifdef DebugLog
341  else LogDebug("HFShower") << "HFShowerLibrary: REJECTED !!!";
342 #endif
343  if (onlyLong && zz >= gpar[4]+gpar[0] && zz <= gpar[4]+gpar[1]) {
344  r1 = G4UniformRand();
345  r2 = G4UniformRand();
346  if (rInside(r) && r1 <= exp(-p*zv) && r2 <= probMax && dfir > gpar[5]){
347  oneHit.position = pos;
348  oneHit.depth = 2;
349  oneHit.time = (tSlice+(pe[i].t())+(fibre->tShift(lpos,2,1)));
350  hit.push_back(oneHit);
351 #ifdef DebugLog
352  LogDebug("HFShower") << "HFShowerLibrary: Final Hit " << nHit
353  <<" position " << (hit[nHit].position) <<" Depth "
354  << (hit[nHit].depth) <<" Time " << tSlice << ":"
355  << pe[i].t() << ":" << fibre->tShift(lpos,2,1)
356  << ":" << (hit[nHit].time);
357 #endif
358  nHit++;
359  }
360  }
361  }
362  }
363 
364 #ifdef DebugLog
365  LogDebug("HFShower") << "HFShowerLibrary: Total Hits " << nHit
366  << " out of " << npe << " PE";
367 #endif
368  if (nHit > npe && !onlyLong)
369  edm::LogWarning("HFShower") << "HFShowerLibrary: Hit buffer " << npe
370  << " smaller than " << nHit << " Hits";
371  return hit;
372 
373 }
374 
375 bool HFShowerLibrary::rInside(double r) {
376 
377  if (r >= rMin && r <= rMax) return true;
378  else return false;
379 }
380 
381 void HFShowerLibrary::getRecord(int type, int record) {
382 
383  int nrc = record-1;
384  int nPhoton = 0;
385  photon.clear();
386  if (type > 0) {
387  hadBranch->SetAddress(&photon);
388  hadBranch->GetEntry(nrc);
389  } else {
390  emBranch->SetAddress(&photon);
391  emBranch->GetEntry(nrc);
392  }
393  nPhoton = photon.size();
394 #ifdef DebugLog
395  LogDebug("HFShower") << "HFShowerLibrary::getRecord: Record " << record
396  << " of type " << type << " with " << nPhoton
397  << " photons";
398  for (int j = 0; j < nPhoton; j++)
399  LogDebug("HFShower") << "Photon " << j << " " << photon[j];
400 #endif
401 }
402 
403 void HFShowerLibrary::loadEventInfo(TBranch* branch) {
404 
405  std::vector<HFShowerLibraryEventInfo> eventInfoCollection;
406  branch->SetAddress(&eventInfoCollection);
407  branch->GetEntry(0);
408  edm::LogInfo("HFShower") << "HFShowerLibrary::loadEventInfo loads "
409  << " EventInfo Collection of size "
410  << eventInfoCollection.size() << " records";
411  totEvents = eventInfoCollection[0].totalEvents();
412  nMomBin = eventInfoCollection[0].numberOfBins();
413  evtPerBin = eventInfoCollection[0].eventsPerBin();
414  libVers = eventInfoCollection[0].showerLibraryVersion();
415  listVersion = eventInfoCollection[0].physListVersion();
416  pmom = eventInfoCollection[0].energyBins();
417  for (unsigned int i=0; i<pmom.size(); i++)
418  pmom[i] *= GeV;
419 }
420 
421 void HFShowerLibrary::interpolate(int type, double pin) {
422 
423 #ifdef DebugLog
424  LogDebug("HFShower") << "HFShowerLibrary:: Interpolate for Energy " <<pin/GeV
425  << " GeV with " << nMomBin << " momentum bins and "
426  << evtPerBin << " entries/bin -- total " << totEvents;
427 #endif
428  int irc[2];
429  double w = 0.;
430  double r = G4UniformRand();
431 
432  if (pin<pmom[0]) {
433  w = pin/pmom[0];
434  irc[1] = int(evtPerBin*r) + 1;
435  irc[0] = 0;
436  } else {
437  for (int j=0; j<nMomBin-1; j++) {
438  if (pin >= pmom[j] && pin < pmom[j+1]) {
439  w = (pin-pmom[j])/(pmom[j+1]-pmom[j]);
440  if (j == nMomBin-2) {
441  irc[1] = int(evtPerBin*0.5*r);
442  } else {
443  irc[1] = int(evtPerBin*r);
444  }
445  irc[1] += (j+1)*evtPerBin + 1;
446  r = G4UniformRand();
447  irc[0] = int(evtPerBin*r) + 1 + j*evtPerBin;
448  if (irc[0]<0) {
449  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[0] = "
450  << irc[0] << " now set to 0";
451  irc[0] = 0;
452  } else if (irc[0] > totEvents) {
453  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[0] = "
454  << irc[0] << " now set to "<< totEvents;
455  irc[0] = totEvents;
456  }
457  }
458  }
459  }
460  if (irc[1]<1) {
461  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[1] = "
462  << irc[1] << " now set to 1";
463  irc[1] = 1;
464  } else if (irc[1] > totEvents) {
465  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[1] = "
466  << irc[1] << " now set to "<< totEvents;
467  irc[1] = totEvents;
468  }
469 
470 #ifdef DebugLog
471  LogDebug("HFShower") << "HFShowerLibrary:: Select records " << irc[0]
472  << " and " << irc[1] << " with weights " << 1-w
473  << " and " << w;
474 #endif
475  pe.clear();
476  npe = 0;
477  int npold = 0;
478  for (int ir=0; ir < 2; ir++) {
479  if (irc[ir]>0) {
480  getRecord (type, irc[ir]);
481  int nPhoton = photon.size();
482  npold += nPhoton;
483  for (int j=0; j<nPhoton; j++) {
484  r = G4UniformRand();
485  if ((ir==0 && r > w) || (ir > 0 && r < w)) {
486  storePhoton (j);
487  }
488  }
489  }
490  }
491 
492  if (npe > npold || (npold == 0 && irc[0] > 0))
493  edm::LogWarning("HFShower") << "HFShowerLibrary: Interpolation Warning =="
494  << " records " << irc[0] << " and " << irc[1]
495  << " gives a buffer of " << npold
496  << " photons and fills " << npe << " *****";
497 #ifdef DebugLog
498  else
499  LogDebug("HFShower") << "HFShowerLibrary: Interpolation == records "
500  << irc[0] << " and " << irc[1] << " gives a "
501  << "buffer of " << npold << " photons and fills "
502  << npe << " PE";
503  for (int j=0; j<npe; j++)
504  LogDebug("HFShower") << "Photon " << j << " " << pe[j];
505 #endif
506 }
507 
508 void HFShowerLibrary::extrapolate(int type, double pin) {
509 
510  int nrec = int(pin/pmom[nMomBin-1]);
511  double w = (pin - pmom[nMomBin-1]*nrec)/pmom[nMomBin-1];
512  nrec++;
513 #ifdef DebugLog
514  LogDebug("HFShower") << "HFShowerLibrary:: Extrapolate for Energy " << pin
515  << " GeV with " << nMomBin << " momentum bins and "
516  << evtPerBin << " entries/bin -- total " << totEvents
517  << " using " << nrec << " records";
518 #endif
519  std::vector<int> irc(nrec);
520 
521  for (int ir=0; ir<nrec; ir++) {
522  double r = G4UniformRand();
523  irc[ir] = int(evtPerBin*0.5*r) +(nMomBin-1)*evtPerBin + 1;
524  if (irc[ir]<1) {
525  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[" << ir
526  << "] = " << irc[ir] << " now set to 1";
527  irc[ir] = 1;
528  } else if (irc[ir] > totEvents) {
529  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[" << ir
530  << "] = " << irc[ir] << " now set to "
531  << totEvents;
532  irc[ir] = totEvents;
533 #ifdef DebugLog
534  } else {
535  LogDebug("HFShower") << "HFShowerLibrary::Extrapolation use irc["
536  << ir << "] = " << irc[ir];
537 #endif
538  }
539  }
540 
541  pe.clear();
542  npe = 0;
543  int npold = 0;
544  for (int ir=0; ir<nrec; ir++) {
545  if (irc[ir]>0) {
546  getRecord (type, irc[ir]);
547  int nPhoton = photon.size();
548  npold += nPhoton;
549  for (int j=0; j<nPhoton; j++) {
550  double r = G4UniformRand();
551  if (ir != nrec-1 || r < w) {
552  storePhoton (j);
553  }
554  }
555 #ifdef DebugLog
556  LogDebug("HFShower") << "HFShowerLibrary: Record [" << ir << "] = "
557  << irc[ir] << " npold = " << npold;
558 #endif
559  }
560  }
561 #ifdef DebugLog
562  edm::LogInfo("HFShower") << "HFShowerLibrary:: uses " << npold << " photons";
563 #endif
564 
565  if (npe > npold || npold == 0)
566  edm::LogWarning("HFShower") << "HFShowerLibrary: Extrapolation Warning == "
567  << nrec << " records " << irc[0] << ", "
568  << irc[1] << ", ... gives a buffer of " <<npold
569  << " photons and fills " << npe
570  << " *****";
571 #ifdef DebugLog
572  else
573  LogDebug("HFShower") << "HFShowerLibrary: Extrapolation == " << nrec
574  << " records " << irc[0] << ", " << irc[1]
575  << ", ... gives a buffer of " << npold
576  << " photons and fills " << npe << " PE";
577  for (int j=0; j<npe; j++)
578  LogDebug("HFShower") << "Photon " << j << " " << pe[j];
579 #endif
580 }
581 
582 void HFShowerLibrary::storePhoton(int j) {
583 
584  pe.push_back(photon[j]);
585 #ifdef DebugLog
586  LogDebug("HFShower") << "HFShowerLibrary: storePhoton " << j << " npe "
587  << npe << " " << pe[npe];
588 #endif
589  npe++;
590 }
591 
592 std::vector<double> HFShowerLibrary::getDDDArray(const std::string & str,
593  const DDsvalues_type & sv,
594  int & nmin) {
595 
596 #ifdef DebugLog
597  LogDebug("HFShower") << "HFShowerLibrary:getDDDArray called for " << str
598  << " with nMin " << nmin;
599 #endif
600  DDValue value(str);
601  if (DDfetch(&sv,value)) {
602 #ifdef DebugLog
603  LogDebug("HFShower") << value;
604 #endif
605  const std::vector<double> & fvec = value.doubles();
606  int nval = fvec.size();
607  if (nmin > 0) {
608  if (nval < nmin) {
609  edm::LogError("HFShower") << "HFShowerLibrary : # of " << str
610  << " bins " << nval << " < " << nmin
611  << " ==> illegal";
612  throw cms::Exception("Unknown", "HFShowerLibrary")
613  << "nval < nmin for array " << str << "\n";
614  }
615  } else {
616  if (nval < 2) {
617  edm::LogError("HFShower") << "HFShowerLibrary : # of " << str
618  << " bins " << nval << " < 2 ==> illegal"
619  << " (nmin=" << nmin << ")";
620  throw cms::Exception("Unknown", "HFShowerLibrary")
621  << "nval < 2 for array " << str << "\n";
622  }
623  }
624  nmin = nval;
625 
626  return fvec;
627  } else {
628  edm::LogError("HFShower") << "HFShowerLibrary : cannot get array " << str;
629  throw cms::Exception("Unknown", "HFShowerLibrary")
630  << "cannot get array " << str << "\n";
631  }
632 }
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:516
type
type
Definition: HCALResponse.h:22
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
i
int i
Definition: DBlmapReader.cc:9
DDValue::doubles
const std::vector< double > & doubles() const
a reference to the double-valued values stored in the given instance of DDValue
Definition: DDValue.cc:151
HFShowerLibrary::getHits
std::vector< Hit > getHits(G4Step *aStep, bool &ok, bool onlyLong=false)
Definition: HFShowerLibrary.cc:187
DDFilteredView::addFilter
void addFilter(const DDFilter &, log_op op=AND)
Definition: DDFilteredView.cc:24
HFShowerLibrary::getRecord
void getRecord(int, int)
Definition: HFShowerLibrary.cc:381
record
JetCorrectorParameters::Record record
Definition: classes.h:11
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
HFFibre::attLength
double attLength(double lambda)
Definition: HFFibre.cc:113
abs
#define abs(x)
Definition: mlp_lapack.h:159
funct::exp
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
DDCompactView
type of data representation of DDCompactView
Definition: DDCompactView.h:81
DDfetch
bool DDfetch(const DDsvalues_type *, DDValue &)
helper for retrieving DDValues from DDsvalues_type *.
Definition: DDsvalues.cc:102
HFShowerLibrary::HFShowerLibrary
HFShowerLibrary(std::string &name, const DDCompactView &cpv, edm::ParameterSet const &p)
Definition: HFShowerLibrary.cc:23
detailsBasic3DVector::z
double double double z
Definition: newBasic3DVector.h:17
HFFibre::zShift
double zShift(G4ThreeVector point, int depth, int fromEndAbs=0)
Definition: HFFibre.cc:143
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
HFFibre::tShift
double tShift(G4ThreeVector point, int depth, int fromEndAbs=0)
Definition: HFFibre.cc:131
DDsvalues_type
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
HFShowerLibrary::loadEventInfo
void loadEventInfo(TBranch *)
Definition: HFShowerLibrary.cc:403
j
int j
Definition: DBlmapReader.cc:9
HFShowerLibrary::gpar
std::vector< double > gpar
Definition: HFShowerLibrary.h:73
HFShowerLibrary::getDDDArray
std::vector< double > getDDDArray(const std::string &, const DDsvalues_type &, int &)
Definition: HFShowerLibrary.cc:592
event
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
HFShowerLibrary::photon
std::vector< HFShowerPhoton > photon
Definition: HFShowerLibrary.h:81
HFShowerLibrary::emBranch
TBranch * emBranch
Definition: HFShowerLibrary.h:64
HFShowerLibrary::pe
std::vector< HFShowerPhoton > pe
Definition: HFShowerLibrary.h:80
HFShowerLibrary::pmom
std::vector< double > pmom
Definition: HFShowerLibrary.h:69
HFShowerLibrary::Hit::position
G4ThreeVector position
Definition: HFShowerLibrary.h:41
align_tpl.filter
tuple filter
USE THIS FOR SKIMMED TRACKS process.p = cms.Path(process.hltLevel1GTSeed*process.skimming*process.offlineBeamSpot*process.TrackRefitter2) OTHERWISE USE THIS.
Definition: align_tpl.py:86
L1TEmulatorMonitor_cff.p
tuple p
Definition: L1TEmulatorMonitor_cff.py:48
DDFilteredView::mergedSpecifics
DDsvalues_type mergedSpecifics() const
Definition: DDFilteredView.cc:66
HFShowerLibrary::initRun
void initRun(G4ParticleTable *theParticleTable)
Definition: HFShowerLibrary.cc:159
HFShowerLibrary::extrapolate
void extrapolate(int, double)
Definition: HFShowerLibrary.cc:508
HFShowerLibrary::storePhoton
void storePhoton(int j)
Definition: HFShowerLibrary.cc:582
DDFilteredView::firstChild
bool firstChild()
set the current node to the first child ...
Definition: DDFilteredView.cc:131
DDSpecificsFilter::setCriteria
void setCriteria(const DDValue &nameVal, comp_op, log_op l=AND, bool asString=true, bool merged=true)
Definition: DDFilter.cc:285
edm::FileInPath::fullPath
std::string fullPath() const
Definition: FileInPath.cc:170
HFShowerLibrary::rInside
bool rInside(double r)
Definition: HFShowerLibrary.cc:375
HFShowerLibrary::interpolate
void interpolate(int, double)
Definition: HFShowerLibrary.cc:421
HFShowerLibrary::hadBranch
TBranch * hadBranch
Definition: HFShowerLibrary.h:64
DDSpecificsFilter
The DDGenericFilter is a runtime-parametrized Filter looking on DDSpecifcs.
Definition: DDFilter.h:37
Generated for CMSSW Reference Manual by doxygen 1.8.5