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