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