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