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