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  photon.clear();
385  if (type > 0) {
386  hadBranch->SetAddress(&photon);
387  hadBranch->GetEntry(nrc);
388  } else {
389  emBranch->SetAddress(&photon);
390  emBranch->GetEntry(nrc);
391  }
392 #ifdef DebugLog
393  int nPhoton = photon.size();
394  LogDebug("HFShower") << "HFShowerLibrary::getRecord: Record " << record
395  << " of type " << type << " with " << nPhoton
396  << " photons";
397  for (int j = 0; j < nPhoton; j++)
398  LogDebug("HFShower") << "Photon " << j << " " << photon[j];
399 #endif
400 }
401 
402 void HFShowerLibrary::loadEventInfo(TBranch* branch) {
403 
404  std::vector<HFShowerLibraryEventInfo> eventInfoCollection;
405  branch->SetAddress(&eventInfoCollection);
406  branch->GetEntry(0);
407  edm::LogInfo("HFShower") << "HFShowerLibrary::loadEventInfo loads "
408  << " EventInfo Collection of size "
409  << eventInfoCollection.size() << " records";
410  totEvents = eventInfoCollection[0].totalEvents();
411  nMomBin = eventInfoCollection[0].numberOfBins();
412  evtPerBin = eventInfoCollection[0].eventsPerBin();
413  libVers = eventInfoCollection[0].showerLibraryVersion();
414  listVersion = eventInfoCollection[0].physListVersion();
415  pmom = eventInfoCollection[0].energyBins();
416  for (unsigned int i=0; i<pmom.size(); i++)
417  pmom[i] *= GeV;
418 }
419 
420 void HFShowerLibrary::interpolate(int type, double pin) {
421 
422 #ifdef DebugLog
423  LogDebug("HFShower") << "HFShowerLibrary:: Interpolate for Energy " <<pin/GeV
424  << " GeV with " << nMomBin << " momentum bins and "
425  << evtPerBin << " entries/bin -- total " << totEvents;
426 #endif
427  int irc[2];
428  double w = 0.;
429  double r = G4UniformRand();
430 
431  if (pin<pmom[0]) {
432  w = pin/pmom[0];
433  irc[1] = int(evtPerBin*r) + 1;
434  irc[0] = 0;
435  } else {
436  for (int j=0; j<nMomBin-1; j++) {
437  if (pin >= pmom[j] && pin < pmom[j+1]) {
438  w = (pin-pmom[j])/(pmom[j+1]-pmom[j]);
439  if (j == nMomBin-2) {
440  irc[1] = int(evtPerBin*0.5*r);
441  } else {
442  irc[1] = int(evtPerBin*r);
443  }
444  irc[1] += (j+1)*evtPerBin + 1;
445  r = G4UniformRand();
446  irc[0] = int(evtPerBin*r) + 1 + j*evtPerBin;
447  if (irc[0]<0) {
448  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[0] = "
449  << irc[0] << " now set to 0";
450  irc[0] = 0;
451  } else if (irc[0] > totEvents) {
452  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[0] = "
453  << irc[0] << " now set to "<< totEvents;
454  irc[0] = totEvents;
455  }
456  }
457  }
458  }
459  if (irc[1]<1) {
460  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[1] = "
461  << irc[1] << " now set to 1";
462  irc[1] = 1;
463  } else if (irc[1] > totEvents) {
464  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[1] = "
465  << irc[1] << " now set to "<< totEvents;
466  irc[1] = totEvents;
467  }
468 
469 #ifdef DebugLog
470  LogDebug("HFShower") << "HFShowerLibrary:: Select records " << irc[0]
471  << " and " << irc[1] << " with weights " << 1-w
472  << " and " << w;
473 #endif
474  pe.clear();
475  npe = 0;
476  int npold = 0;
477  for (int ir=0; ir < 2; ir++) {
478  if (irc[ir]>0) {
479  getRecord (type, irc[ir]);
480  int nPhoton = photon.size();
481  npold += nPhoton;
482  for (int j=0; j<nPhoton; j++) {
483  r = G4UniformRand();
484  if ((ir==0 && r > w) || (ir > 0 && r < w)) {
485  storePhoton (j);
486  }
487  }
488  }
489  }
490 
491  if (npe > npold || (npold == 0 && irc[0] > 0))
492  edm::LogWarning("HFShower") << "HFShowerLibrary: Interpolation Warning =="
493  << " records " << irc[0] << " and " << irc[1]
494  << " gives a buffer of " << npold
495  << " photons and fills " << npe << " *****";
496 #ifdef DebugLog
497  else
498  LogDebug("HFShower") << "HFShowerLibrary: Interpolation == records "
499  << irc[0] << " and " << irc[1] << " gives a "
500  << "buffer of " << npold << " photons and fills "
501  << npe << " PE";
502  for (int j=0; j<npe; j++)
503  LogDebug("HFShower") << "Photon " << j << " " << pe[j];
504 #endif
505 }
506 
507 void HFShowerLibrary::extrapolate(int type, double pin) {
508 
509  int nrec = int(pin/pmom[nMomBin-1]);
510  double w = (pin - pmom[nMomBin-1]*nrec)/pmom[nMomBin-1];
511  nrec++;
512 #ifdef DebugLog
513  LogDebug("HFShower") << "HFShowerLibrary:: Extrapolate for Energy " << pin
514  << " GeV with " << nMomBin << " momentum bins and "
515  << evtPerBin << " entries/bin -- total " << totEvents
516  << " using " << nrec << " records";
517 #endif
518  std::vector<int> irc(nrec);
519 
520  for (int ir=0; ir<nrec; ir++) {
521  double r = G4UniformRand();
522  irc[ir] = int(evtPerBin*0.5*r) +(nMomBin-1)*evtPerBin + 1;
523  if (irc[ir]<1) {
524  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[" << ir
525  << "] = " << irc[ir] << " now set to 1";
526  irc[ir] = 1;
527  } else if (irc[ir] > totEvents) {
528  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[" << ir
529  << "] = " << irc[ir] << " now set to "
530  << totEvents;
531  irc[ir] = totEvents;
532 #ifdef DebugLog
533  } else {
534  LogDebug("HFShower") << "HFShowerLibrary::Extrapolation use irc["
535  << ir << "] = " << irc[ir];
536 #endif
537  }
538  }
539 
540  pe.clear();
541  npe = 0;
542  int npold = 0;
543  for (int ir=0; ir<nrec; ir++) {
544  if (irc[ir]>0) {
545  getRecord (type, irc[ir]);
546  int nPhoton = photon.size();
547  npold += nPhoton;
548  for (int j=0; j<nPhoton; j++) {
549  double r = G4UniformRand();
550  if (ir != nrec-1 || r < w) {
551  storePhoton (j);
552  }
553  }
554 #ifdef DebugLog
555  LogDebug("HFShower") << "HFShowerLibrary: Record [" << ir << "] = "
556  << irc[ir] << " npold = " << npold;
557 #endif
558  }
559  }
560 #ifdef DebugLog
561  edm::LogInfo("HFShower") << "HFShowerLibrary:: uses " << npold << " photons";
562 #endif
563 
564  if (npe > npold || npold == 0)
565  edm::LogWarning("HFShower") << "HFShowerLibrary: Extrapolation Warning == "
566  << nrec << " records " << irc[0] << ", "
567  << irc[1] << ", ... gives a buffer of " <<npold
568  << " photons and fills " << npe
569  << " *****";
570 #ifdef DebugLog
571  else
572  LogDebug("HFShower") << "HFShowerLibrary: Extrapolation == " << nrec
573  << " records " << irc[0] << ", " << irc[1]
574  << ", ... gives a buffer of " << npold
575  << " photons and fills " << npe << " PE";
576  for (int j=0; j<npe; j++)
577  LogDebug("HFShower") << "Photon " << j << " " << pe[j];
578 #endif
579 }
580 
581 void HFShowerLibrary::storePhoton(int j) {
582 
583  pe.push_back(photon[j]);
584 #ifdef DebugLog
585  LogDebug("HFShower") << "HFShowerLibrary: storePhoton " << j << " npe "
586  << npe << " " << pe[npe];
587 #endif
588  npe++;
589 }
590 
591 std::vector<double> HFShowerLibrary::getDDDArray(const std::string & str,
592  const DDsvalues_type & sv,
593  int & nmin) {
594 
595 #ifdef DebugLog
596  LogDebug("HFShower") << "HFShowerLibrary:getDDDArray called for " << str
597  << " with nMin " << nmin;
598 #endif
599  DDValue value(str);
600  if (DDfetch(&sv,value)) {
601 #ifdef DebugLog
602  LogDebug("HFShower") << value;
603 #endif
604  const std::vector<double> & fvec = value.doubles();
605  int nval = fvec.size();
606  if (nmin > 0) {
607  if (nval < nmin) {
608  edm::LogError("HFShower") << "HFShowerLibrary : # of " << str
609  << " bins " << nval << " < " << nmin
610  << " ==> illegal";
611  throw cms::Exception("Unknown", "HFShowerLibrary")
612  << "nval < nmin for array " << str << "\n";
613  }
614  } else {
615  if (nval < 2) {
616  edm::LogError("HFShower") << "HFShowerLibrary : # of " << str
617  << " bins " << nval << " < 2 ==> illegal"
618  << " (nmin=" << nmin << ")";
619  throw cms::Exception("Unknown", "HFShowerLibrary")
620  << "nval < 2 for array " << str << "\n";
621  }
622  }
623  nmin = nval;
624 
625  return fvec;
626  } else {
627  edm::LogError("HFShower") << "HFShowerLibrary : cannot get array " << str;
628  throw cms::Exception("Unknown", "HFShowerLibrary")
629  << "cannot get array " << str << "\n";
630  }
631 }
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:501
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:134
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
create_public_lumi_plots.exp
tuple exp
Definition: create_public_lumi_plots.py:1093
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
alcazmumu_cfi.filter
tuple filter
Definition: alcazmumu_cfi.py:34
DDCompactView
type of data representation of DDCompactView
Definition: DDCompactView.h:77
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:402
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:591
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
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:507
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
HFShowerLibrary::storePhoton
void storePhoton(int j)
Definition: HFShowerLibrary.cc:581
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:171
HFShowerLibrary::rInside
bool rInside(double r)
Definition: HFShowerLibrary.cc:375
HFShowerLibrary::interpolate
void interpolate(int, double)
Definition: HFShowerLibrary.cc:420
w
T w() const
Definition: newBasic3DVector.h:234
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