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,DDSpecificsFilter::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  std::vector<HFShowerLibrary::Hit> hit;
215  ok = false;
216  if (parCode == pi0PDG || parCode == etaPDG || parCode == nuePDG ||
217  parCode == numuPDG || parCode == nutauPDG || parCode == anuePDG ||
218  parCode == anumuPDG || parCode == anutauPDG || parCode == geantinoPDG)
219  return hit;
220  ok = true;
221 
222  double tSlice = (postStepPoint->GetGlobalTime())/nanosecond;
223  double pin = preStepPoint->GetTotalEnergy();
224  double pz = momDir.z();
225  double zint = hitPoint.z();
226 
227  // if particle moves from interaction point or "backwards (halo)
228  int backward = 0;
229  if (pz * zint < 0.) backward = 1;
230 
231  double sphi = sin(momDir.phi());
232  double cphi = cos(momDir.phi());
233  double ctheta = cos(momDir.theta());
234  double stheta = sin(momDir.theta());
235 
236 #ifdef DebugLog
237  G4ThreeVector localPos = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
238  double zoff = localPos.z() + 0.5*gpar[1];
239  // if (zoff < 0) zoff = 0;
240  edm::LogInfo("HFShower") << "HFShowerLibrary: getHits " << partType
241  << " of energy " << pin/GeV << " GeV"
242  << " dir.orts " << momDir.x() << ", " <<momDir.y()
243  << ", " << momDir.z() << " Pos x,y,z = "
244  << hitPoint.x() << "," << hitPoint.y() << ","
245  << hitPoint.z() << " (" << zoff
246  << ") sphi,cphi,stheta,ctheta = " << sphi
247  << "," << cphi << ", " << stheta << "," << ctheta;
248 #endif
249 
250  if (parCode == emPDG || parCode == epPDG || parCode == gammaPDG ) {
251  if (pin<pmom[nMomBin-1]) {
252  interpolate(0, pin);
253  } else {
254  extrapolate(0, pin);
255  }
256  } else {
257  if (pin<pmom[nMomBin-1]) {
258  interpolate(1, pin);
259  } else {
260  extrapolate(1, pin);
261  }
262  }
263 
264  int nHit = 0;
265  HFShowerLibrary::Hit oneHit;
266  for (int i = 0; i < npe; i++) {
267  double zv = std::abs(pe[i].z()); // abs local z
268 #ifdef DebugLog
269  edm::LogInfo("HFShower") << "HFShowerLibrary: Hit " << i << " " << pe[i] << " zv " << zv;
270 #endif
271  if (zv <= gpar[1] && pe[i].lambda() > 0 &&
272  (pe[i].z() >= 0 || (zv > gpar[0] && (!onlyLong)))) {
273  int depth = 1;
274  if (onlyLong) {
275  } else if (backward == 0) { // fully valid only for "front" particles
276  if (pe[i].z() < 0) depth = 2;// with "front"-simulated shower lib.
277  } else { // for "backward" particles - almost equal
278  double r = G4UniformRand(); // share between L and S fibers
279  if (r > 0.5) depth = 2;
280  }
281 
282 
283  // Updated coordinate transformation from local
284  // back to global using two Euler angles: phi and theta
285  double pex = pe[i].x();
286  double pey = pe[i].y();
287 
288  double xx = pex*ctheta*cphi - pey*sphi + zv*stheta*cphi;
289  double yy = pex*ctheta*sphi + pey*cphi + zv*stheta*sphi;
290  double zz = -pex*stheta + zv*ctheta;
291 
292  G4ThreeVector pos = hitPoint + G4ThreeVector(xx,yy,zz);
293  zv = std::abs(pos.z()) - gpar[4] - 0.5*gpar[1];
294  G4ThreeVector lpos = G4ThreeVector(pos.x(),pos.y(),zv);
295 
296  zv = fibre->zShift(lpos,depth,0); // distance to PMT !
297 
298  double r = pos.perp();
299  double p = fibre->attLength(pe[i].lambda());
300  double fi = pos.phi();
301  if (fi < 0) fi += twopi;
302  int isect = int(fi/dphi) + 1;
303  isect = (isect + 1) / 2;
304  double dfi = ((isect*2-1)*dphi - fi);
305  if (dfi < 0) dfi = -dfi;
306  double dfir = r * sin(dfi);
307 #ifdef DebugLog
308  edm::LogInfo("HFShower") << "HFShowerLibrary: Position shift " << xx
309  << ", " << yy << ", " << zz << ": " << pos
310  << " R " << r << " Phi " << fi << " Section "
311  << isect << " R*Dfi " << dfir << " Dist " << zv;
312 #endif
313  zz = std::abs(pos.z());
314  double r1 = G4UniformRand();
315  double r2 = G4UniformRand();
316  double r3 = -9999.;
317  if (backward) r3 = G4UniformRand();
318  if (!applyFidCut) dfir += gpar[5];
319 
320 #ifdef DebugLog
321  edm::LogInfo("HFShower") << "HFShowerLibrary: rLimits " << rInside(r)
322  << " attenuation " << r1 <<":" << exp(-p*zv)
323  << " r2 " << r2 << " r3 " << r3 << " rDfi "
324  << gpar[5] << " zz "
325  << zz << " zLim " << gpar[4] << ":"
326  << gpar[4]+gpar[1] << "\n"
327  << " rInside(r) :" << rInside(r)
328  << " r1 <= exp(-p*zv) :" << (r1 <= exp(-p*zv))
329  << " r2 <= probMax :" << (r2 <= probMax*weight)
330  << " r3 <= backProb :" << (r3 <= backProb)
331  << " dfir > gpar[5] :" << (dfir > gpar[5])
332  << " zz >= gpar[4] :" << (zz >= gpar[4])
333  << " zz <= gpar[4]+gpar[1] :"
334  << (zz <= gpar[4]+gpar[1]);
335 #endif
336  if (rInside(r) && r1 <= exp(-p*zv) && r2 <= probMax*weight &&
337  dfir > gpar[5] && zz >= gpar[4] && zz <= gpar[4]+gpar[1] &&
338  r3 <= backProb && (depth != 2 || zz >= gpar[4]+gpar[0])) {
339  oneHit.position = pos;
340  oneHit.depth = depth;
341  oneHit.time = (tSlice+(pe[i].t())+(fibre->tShift(lpos,depth,1)));
342  hit.push_back(oneHit);
343 #ifdef DebugLog
344  edm::LogInfo("HFShower") << "HFShowerLibrary: Final Hit " << nHit
345  <<" position " << (hit[nHit].position)
346  << " Depth " << (hit[nHit].depth) <<" Time "
347  << tSlice << ":" << pe[i].t() << ":"
348  << fibre->tShift(lpos,depth,1) << ":"
349  << (hit[nHit].time);
350 #endif
351  nHit++;
352  }
353 #ifdef DebugLog
354  else LogDebug("HFShower") << "HFShowerLibrary: REJECTED !!!";
355 #endif
356  if (onlyLong && zz >= gpar[4]+gpar[0] && zz <= gpar[4]+gpar[1]) {
357  r1 = G4UniformRand();
358  r2 = G4UniformRand();
359  if (rInside(r) && r1 <= exp(-p*zv) && r2 <= probMax && dfir > gpar[5]){
360  oneHit.position = pos;
361  oneHit.depth = 2;
362  oneHit.time = (tSlice+(pe[i].t())+(fibre->tShift(lpos,2,1)));
363  hit.push_back(oneHit);
364 #ifdef DebugLog
365  edm::LogInfo("HFShower") << "HFShowerLibrary: Final Hit " << nHit
366  << " position " << (hit[nHit].position)
367  << " Depth " << (hit[nHit].depth) <<" Time "
368  << tSlice << ":" << pe[i].t() << ":"
369  << fibre->tShift(lpos,2,1) << ":"
370  << (hit[nHit].time);
371 #endif
372  nHit++;
373  }
374  }
375  }
376  }
377 
378 #ifdef DebugLog
379  edm::LogInfo("HFShower") << "HFShowerLibrary: Total Hits " << nHit
380  << " out of " << npe << " PE";
381 #endif
382  if (nHit > npe && !onlyLong)
383  edm::LogWarning("HFShower") << "HFShowerLibrary: Hit buffer " << npe
384  << " smaller than " << nHit << " Hits";
385  return hit;
386 
387 }
388 
389 bool HFShowerLibrary::rInside(double r) {
390 
391  if (r >= rMin && r <= rMax) return true;
392  else return false;
393 }
394 
395 void HFShowerLibrary::getRecord(int type, int record) {
396 
397  int nrc = record-1;
398  photon.clear();
399  photo->clear();
400  if (type > 0) {
401  if (newForm) {
402  hadBranch->SetAddress(&photo);
403  hadBranch->GetEntry(nrc+totEvents);
404  } else {
405  hadBranch->SetAddress(&photon);
406  hadBranch->GetEntry(nrc);
407  }
408  } else {
409  if (newForm) {
410  emBranch->SetAddress(&photo);
411  } else {
412  emBranch->SetAddress(&photon);
413  }
414  emBranch->GetEntry(nrc);
415  }
416 #ifdef DebugLog
417  int nPhoton = (newForm) ? photo->size() : photon.size();
418  LogDebug("HFShower") << "HFShowerLibrary::getRecord: Record " << record
419  << " of type " << type << " with " << nPhoton
420  << " photons";
421  for (int j = 0; j < nPhoton; j++)
422  if (newForm) LogDebug("HFShower") << "Photon " << j << " " << photo->at(j);
423  else LogDebug("HFShower") << "Photon " << j << " " << photon[j];
424 #endif
425 }
426 
427 void HFShowerLibrary::loadEventInfo(TBranch* branch) {
428 
429  if (branch) {
430  std::vector<HFShowerLibraryEventInfo> eventInfoCollection;
431  branch->SetAddress(&eventInfoCollection);
432  branch->GetEntry(0);
433  edm::LogInfo("HFShower") << "HFShowerLibrary::loadEventInfo loads "
434  << " EventInfo Collection of size "
435  << eventInfoCollection.size() << " records";
436  totEvents = eventInfoCollection[0].totalEvents();
437  nMomBin = eventInfoCollection[0].numberOfBins();
438  evtPerBin = eventInfoCollection[0].eventsPerBin();
439  libVers = eventInfoCollection[0].showerLibraryVersion();
440  listVersion = eventInfoCollection[0].physListVersion();
441  pmom = eventInfoCollection[0].energyBins();
442  } else {
443  edm::LogInfo("HFShower") << "HFShowerLibrary::loadEventInfo loads "
444  << " EventInfo from hardwired numbers";
445  nMomBin = 16;
446  evtPerBin = 5000;
447  totEvents = nMomBin*evtPerBin;
448  libVers = 1.1;
449  listVersion = 3.6;
450  pmom = {2,3,5,7,10,15,20,30,50,75,100,150,250,350,500,1000};
451  }
452  for (int i=0; i<nMomBin; i++)
453  pmom[i] *= GeV;
454 }
455 
456 void HFShowerLibrary::interpolate(int type, double pin) {
457 
458 #ifdef DebugLog
459  LogDebug("HFShower") << "HFShowerLibrary:: Interpolate for Energy " <<pin/GeV
460  << " GeV with " << nMomBin << " momentum bins and "
461  << evtPerBin << " entries/bin -- total " << totEvents;
462 #endif
463  int irc[2];
464  double w = 0.;
465  double r = G4UniformRand();
466 
467  if (pin<pmom[0]) {
468  w = pin/pmom[0];
469  irc[1] = int(evtPerBin*r) + 1;
470  irc[0] = 0;
471  } else {
472  for (int j=0; j<nMomBin-1; j++) {
473  if (pin >= pmom[j] && pin < pmom[j+1]) {
474  w = (pin-pmom[j])/(pmom[j+1]-pmom[j]);
475  if (j == nMomBin-2) {
476  irc[1] = int(evtPerBin*0.5*r);
477  } else {
478  irc[1] = int(evtPerBin*r);
479  }
480  irc[1] += (j+1)*evtPerBin + 1;
481  r = G4UniformRand();
482  irc[0] = int(evtPerBin*r) + 1 + j*evtPerBin;
483  if (irc[0]<0) {
484  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[0] = "
485  << irc[0] << " now set to 0";
486  irc[0] = 0;
487  } else if (irc[0] > totEvents) {
488  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[0] = "
489  << irc[0] << " now set to "<< totEvents;
490  irc[0] = totEvents;
491  }
492  }
493  }
494  }
495  if (irc[1]<1) {
496  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[1] = "
497  << irc[1] << " now set to 1";
498  irc[1] = 1;
499  } else if (irc[1] > totEvents) {
500  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[1] = "
501  << irc[1] << " now set to "<< totEvents;
502  irc[1] = totEvents;
503  }
504 
505 #ifdef DebugLog
506  LogDebug("HFShower") << "HFShowerLibrary:: Select records " << irc[0]
507  << " and " << irc[1] << " with weights " << 1-w
508  << " and " << w;
509 #endif
510  pe.clear();
511  npe = 0;
512  int npold = 0;
513  for (int ir=0; ir < 2; ir++) {
514  if (irc[ir]>0) {
515  getRecord (type, irc[ir]);
516  int nPhoton = (newForm) ? photo->size() : photon.size();
517  npold += nPhoton;
518  for (int j=0; j<nPhoton; j++) {
519  r = G4UniformRand();
520  if ((ir==0 && r > w) || (ir > 0 && r < w)) {
521  storePhoton (j);
522  }
523  }
524  }
525  }
526 
527  if ((npe > npold || (npold == 0 && irc[0] > 0)) && !(npe == 0 && npold == 0))
528  edm::LogWarning("HFShower") << "HFShowerLibrary: Interpolation Warning =="
529  << " records " << irc[0] << " and " << irc[1]
530  << " gives a buffer of " << npold
531  << " photons and fills " << npe << " *****";
532 #ifdef DebugLog
533  else
534  LogDebug("HFShower") << "HFShowerLibrary: Interpolation == records "
535  << irc[0] << " and " << irc[1] << " gives a "
536  << "buffer of " << npold << " photons and fills "
537  << npe << " PE";
538  for (int j=0; j<npe; j++)
539  LogDebug("HFShower") << "Photon " << j << " " << pe[j];
540 #endif
541 }
542 
543 void HFShowerLibrary::extrapolate(int type, double pin) {
544 
545  int nrec = int(pin/pmom[nMomBin-1]);
546  double w = (pin - pmom[nMomBin-1]*nrec)/pmom[nMomBin-1];
547  nrec++;
548 #ifdef DebugLog
549  LogDebug("HFShower") << "HFShowerLibrary:: Extrapolate for Energy " << pin
550  << " GeV with " << nMomBin << " momentum bins and "
551  << evtPerBin << " entries/bin -- total " << totEvents
552  << " using " << nrec << " records";
553 #endif
554  std::vector<int> irc(nrec);
555 
556  for (int ir=0; ir<nrec; ir++) {
557  double r = G4UniformRand();
558  irc[ir] = int(evtPerBin*0.5*r) +(nMomBin-1)*evtPerBin + 1;
559  if (irc[ir]<1) {
560  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[" << ir
561  << "] = " << irc[ir] << " now set to 1";
562  irc[ir] = 1;
563  } else if (irc[ir] > totEvents) {
564  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[" << ir
565  << "] = " << irc[ir] << " now set to "
566  << totEvents;
567  irc[ir] = totEvents;
568 #ifdef DebugLog
569  } else {
570  LogDebug("HFShower") << "HFShowerLibrary::Extrapolation use irc["
571  << ir << "] = " << irc[ir];
572 #endif
573  }
574  }
575 
576  pe.clear();
577  npe = 0;
578  int npold = 0;
579  for (int ir=0; ir<nrec; ir++) {
580  if (irc[ir]>0) {
581  getRecord (type, irc[ir]);
582  int nPhoton = (newForm) ? photo->size() : photon.size();
583  npold += nPhoton;
584  for (int j=0; j<nPhoton; j++) {
585  double r = G4UniformRand();
586  if (ir != nrec-1 || r < w) {
587  storePhoton (j);
588  }
589  }
590 #ifdef DebugLog
591  LogDebug("HFShower") << "HFShowerLibrary: Record [" << ir << "] = "
592  << irc[ir] << " npold = " << npold;
593 #endif
594  }
595  }
596 #ifdef DebugLog
597  LogDebug("HFShower") << "HFShowerLibrary:: uses " << npold << " photons";
598 #endif
599 
600  if (npe > npold || npold == 0)
601  edm::LogWarning("HFShower") << "HFShowerLibrary: Extrapolation Warning == "
602  << nrec << " records " << irc[0] << ", "
603  << irc[1] << ", ... gives a buffer of " <<npold
604  << " photons and fills " << npe
605  << " *****";
606 #ifdef DebugLog
607  else
608  LogDebug("HFShower") << "HFShowerLibrary: Extrapolation == " << nrec
609  << " records " << irc[0] << ", " << irc[1]
610  << ", ... gives a buffer of " << npold
611  << " photons and fills " << npe << " PE";
612  for (int j=0; j<npe; j++)
613  LogDebug("HFShower") << "Photon " << j << " " << pe[j];
614 #endif
615 }
616 
617 void HFShowerLibrary::storePhoton(int j) {
618 
619  if (newForm) pe.push_back(photo->at(j));
620  else pe.push_back(photon[j]);
621 #ifdef DebugLog
622  LogDebug("HFShower") << "HFShowerLibrary: storePhoton " << j << " npe "
623  << npe << " " << pe[npe];
624 #endif
625  npe++;
626 }
627 
628 std::vector<double> HFShowerLibrary::getDDDArray(const std::string & str,
629  const DDsvalues_type & sv,
630  int & nmin) {
631 
632 #ifdef DebugLog
633  LogDebug("HFShower") << "HFShowerLibrary:getDDDArray called for " << str
634  << " with nMin " << nmin;
635 #endif
636  DDValue value(str);
637  if (DDfetch(&sv,value)) {
638 #ifdef DebugLog
639  LogDebug("HFShower") << value;
640 #endif
641  const std::vector<double> & fvec = value.doubles();
642  int nval = fvec.size();
643  if (nmin > 0) {
644  if (nval < nmin) {
645  edm::LogError("HFShower") << "HFShowerLibrary : # of " << str
646  << " bins " << nval << " < " << nmin
647  << " ==> illegal";
648  throw cms::Exception("Unknown", "HFShowerLibrary")
649  << "nval < nmin for array " << str << "\n";
650  }
651  } else {
652  if (nval < 2) {
653  edm::LogError("HFShower") << "HFShowerLibrary : # of " << str
654  << " bins " << nval << " < 2 ==> illegal"
655  << " (nmin=" << nmin << ")";
656  throw cms::Exception("Unknown", "HFShowerLibrary")
657  << "nval < 2 for array " << str << "\n";
658  }
659  }
660  nmin = nval;
661 
662  return fvec;
663  } else {
664  edm::LogError("HFShower") << "HFShowerLibrary : cannot get array " << str;
665  throw cms::Exception("Unknown", "HFShowerLibrary")
666  << "cannot get array " << str << "\n";
667  }
668 }
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:139
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:170
GeV
const double GeV
Definition: MathUtil.h:16
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:395
w
const double w
Definition: UKUtility.cc:23
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:81
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:80
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
HFShowerPhotonCollection
std::vector< HFShowerPhoton > HFShowerPhotonCollection
Definition: HFShowerPhoton.h:40
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:427
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:628
HFShowerLibrary::photo
HFShowerPhotonCollection * photo
Definition: HFShowerLibrary.h:82
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:65
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:169
HFShowerLibrary::extrapolate
void extrapolate(int, double)
Definition: HFShowerLibrary.cc:543
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
HFShowerLibrary::storePhoton
void storePhoton(int j)
Definition: HFShowerLibrary.cc:617
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
histoStyle.weight
int weight
Definition: histoStyle.py:50
HFShowerLibrary::photon
HFShowerPhotonCollection photon
Definition: HFShowerLibrary.h:83
HFShowerLibrary::rInside
bool rInside(double r)
Definition: HFShowerLibrary.cc:389
HFShowerLibrary::interpolate
void interpolate(int, double)
Definition: HFShowerLibrary.cc:456
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