CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
FastHFShowerLibrary.cc
Go to the documentation of this file.
1 // File: HFShowerLibrary.cc
3 // Description: Shower library for Very forward hadron calorimeter
5 
6 #include "FastSimulation/ShowerDevelopment/interface/FastHFShowerLibrary.h"
7 #include "FastSimulation/Event/interface/FSimEvent.h"
8 #include "FastSimulation/Event/interface/FSimTrack.h"
9 #include "SimG4CMS/Calo/interface/HFFibreFiducial.h"
10 #include "SimDataFormats/CaloHit/interface/HFShowerLibraryEventInfo.h"
11 #include "DetectorDescription/Core/interface/DDFilter.h"
12 #include "DetectorDescription/Core/interface/DDFilteredView.h"
13 #include "DetectorDescription/Core/interface/DDValue.h"
14 #include "FWCore/Framework/interface/ESTransientHandle.h"
15 #include "Geometry/Records/interface/IdealGeometryRecord.h"
16 #include "FWCore/Utilities/interface/Exception.h"
17 
18 #include "Randomize.hh"
19 #include "CLHEP/Units/GlobalSystemOfUnits.h"
20 #include "CLHEP/Units/GlobalPhysicalConstants.h"
21 
22 // STL headers
23 #include <vector>
24 #include <iostream>
25 
26 //#define DebugLog
27 
28 FastHFShowerLibrary::FastHFShowerLibrary(edm::ParameterSet const & p) : fibre(0),
29  hf(0),
30  emBranch(0),
31  hadBranch(0),
32  numberingScheme(0),
33  numberingFromDDD(0),
34  npe(0),
35  photo(0) {
36 
37  edm::ParameterSet m_HS = p.getParameter<edm::ParameterSet>("HFShowerLibrary");
38  edm::FileInPath fp = m_HS.getParameter<edm::FileInPath>("FileName");
39  std::string pTreeName = fp.fullPath();
40  std::string emName = m_HS.getParameter<std::string>("TreeEMID");
41  std::string hadName = m_HS.getParameter<std::string>("TreeHadID");
42  std::string branchEvInfo = m_HS.getUntrackedParameter<std::string>("BranchEvt","HFShowerLibraryEventInfos_hfshowerlib_HFShowerLibraryEventInfo");
43  std::string branchPre = m_HS.getUntrackedParameter<std::string>("BranchPre","HFShowerPhotons_hfshowerlib_");
44  std::string branchPost = m_HS.getUntrackedParameter<std::string>("BranchPost","_R.obj");
45 
46  probMax = m_HS.getParameter<double>("ProbMax");
47  backProb = m_HS.getParameter<double>("BackProbability");
48  verbose = m_HS.getUntrackedParameter<bool>("Verbosity",false);
49  applyFidCut = m_HS.getParameter<bool>("ApplyFiducialCut");
50  cFibre = c_light*(m_HS.getParameter<double>("CFibre"));
51 
52  if (pTreeName.find(".") == 0) pTreeName.erase(0,2);
53  const char* nTree = pTreeName.c_str();
54  hf = TFile::Open(nTree);
55 
56  if (!hf->IsOpen()) {
57  edm::LogError("FastCalorimetry") << "HFShowerLibrary: opening " << nTree
58  << " failed";
59  throw cms::Exception("Unknown", "HFShowerLibrary")
60  << "Opening of " << pTreeName << " fails\n";
61  } else {
62  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: opening " << nTree
63  << " successfully";
64  }
65 
66  newForm = (branchEvInfo == "");
67  TTree* event(0);
68  if (newForm) event = (TTree *) hf ->Get("HFSimHits");
69  else event = (TTree *) hf ->Get("Events");
70  if (event) {
71  TBranch *evtInfo(0);
72  if (!newForm) {
73  std::string info = branchEvInfo + branchPost;
74  evtInfo = event->GetBranch(info.c_str());
75  }
76  if (evtInfo || newForm) {
77  loadEventInfo(evtInfo);
78  } else {
79  edm::LogError("FastCalorimetry") << "HFShowerLibrary: HFShowerLibrayEventInfo"
80  << " Branch does not exist in Event";
81  throw cms::Exception("Unknown", "HFShowerLibrary")
82  << "Event information absent\n";
83  }
84  } else {
85  edm::LogError("FastCalorimetry") << "HFShowerLibrary: Events Tree does not "
86  << "exist";
87  throw cms::Exception("Unknown", "HFShowerLibrary")
88  << "Events tree absent\n";
89  }
90 
91  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: Library " << libVers
92  << " ListVersion " << listVersion
93  << " Events Total " << totEvents << " and "
94  << evtPerBin << " per bin";
95 
96  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: Energies (GeV) with "
97  << nMomBin << " bins";
98 
99  for (int i=0; i<nMomBin; i++)
100  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: pmom[" << i << "] = "
101  << pmom[i]/GeV << " GeV";
102 
103  std::string nameBr = branchPre + emName + branchPost;
104  emBranch = event->GetBranch(nameBr.c_str());
105  if (verbose) emBranch->Print();
106  nameBr = branchPre + hadName + branchPost;
107  hadBranch = event->GetBranch(nameBr.c_str());
108  if (verbose) hadBranch->Print();
109  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary:Branch " << emName
110  << " has " << emBranch->GetEntries()
111  << " entries and Branch " << hadName
112  << " has " << hadBranch->GetEntries()
113  << " entries";
114 
115  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary::No packing information -"
116  << " Assume x, y, z are not in packed form";
117 
118  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: Maximum probability cut off "
119  << probMax << " Back propagation of light prob. "
120  << backProb ;
121 }
122 
123 FastHFShowerLibrary::~FastHFShowerLibrary() {
124  if (hf) hf->Close();
125  if (fibre) delete fibre;
126  if (photo) delete photo;
127  if (numberingFromDDD) delete numberingFromDDD;
128  if (numberingScheme) delete numberingScheme;
129 }
130 
131 void const FastHFShowerLibrary::initHFShowerLibrary(const edm::EventSetup& iSetup) {
132 
133  edm::LogInfo("FastCalorimetry") << "initHFShowerLibrary::initialization";
134 
135  edm::ESTransientHandle<DDCompactView> cpv;
136  iSetup.get<IdealGeometryRecord>().get(cpv);
137 
138  std::string name = "HcalHits";
139  std::string attribute = "ReadOutName";
140  std::string value = name;
141  DDSpecificsFilter filter;
142  DDValue ddv(attribute,value,0);
143  filter.setCriteria(ddv,DDSpecificsFilter::equals);
144  DDFilteredView fv(*cpv);
145  fv.addFilter(filter);
146  bool dodet = fv.firstChild();
147  if (dodet) {
148  DDsvalues_type sv(fv.mergedSpecifics());
149  // Radius (minimum and maximum)
150  int nR = -1;
151  std::vector<double> rTable = getDDDArray("rTable",sv,nR);
152  rMin = rTable[0];
153  rMax = rTable[nR-1];
154  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: rMIN " << rMin/cm
155  << " cm and rMax " << rMax/cm;
156  // Delta phi
157  int nEta = -1;
158  std::vector<double> etaTable = getDDDArray("etaTable",sv,nEta);
159  int nPhi = nEta + nR - 2;
160  std::vector<double> phibin = getDDDArray("phibin",sv,nPhi);
161  dphi = phibin[nEta-1];
162  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: (Half) Phi Width of wedge "
163  << dphi/deg;
164 
165  // Special Geometry parameters
166  int ngpar = 7;
167  gpar = getDDDArray("gparHF",sv,ngpar);
168  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: " << ngpar << " gpar (cm)";
169 
170  for (int ig=0; ig<ngpar; ig++)
171  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: gpar[" << ig << "] = "
172  << gpar[ig]/cm << " cm";
173  } else {
174  edm::LogError("FastCalorimetry") << "HFShowerLibrary: cannot get filtered "
175  << " view for " << attribute << " matching "
176  << name;
177  throw cms::Exception("Unknown", "HFShowerLibrary")
178  << "cannot match " << attribute << " to " << name <<"\n";
179  }
180 
181  initRun();
182  fibre = new FastHFFibre(name, *cpv, cFibre);
183  photo = new HFShowerPhotonCollection;
184  numberingFromDDD = new HcalNumberingFromDDD(name, *cpv);
185  numberingScheme = new HcalNumberingScheme();
186 }
187 
188 void FastHFShowerLibrary::initRun() {
189 
190  geantinoPDG = 0; gammaPDG = 22;
191  emPDG = 11; epPDG = -11; nuePDG = 12; anuePDG = -12;
192  numuPDG = 14; anumuPDG = -14; nutauPDG = 16; anutauPDG = -16;
193  pi0PDG = 111; etaPDG = 221;
194 
195 #ifdef DebugLog
196  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: Particle codes for e- = "
197  << emPDG << ", e+ = " << epPDG << ", gamma = "
198  << gammaPDG << ", pi0 = " << pi0PDG << ", eta = "
199  << etaPDG << ", geantino = " << geantinoPDG
200  << "\n nu_e = " << nuePDG << ", nu_mu = "
201  << numuPDG << ", nu_tau = " << nutauPDG
202  << ", anti_nu_e = " << anuePDG << ", anti_nu_mu = "
203  << anumuPDG << ", anti_nu_tau = " << anutauPDG;
204 #endif
205 }
206 
207 void FastHFShowerLibrary::recoHFShowerLibrary(const FSimTrack& myTrack) {
208 
209 #ifdef DebugLog
210  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: recoHFShowerLibrary ";
211 #endif
212 
213  if(!myTrack.onVFcal()) {
214 #ifdef DebugLog
215  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: we should not be here ";
216 #endif
217  }
218 
219  hitMap.clear();
220  double eGen = 1000.*myTrack.vfcalEntrance().e(); // energy in [MeV]
221  double delZv = (myTrack.vfcalEntrance().vertex().Z()>0.0) ? 50.0 : -50.0;
222  G4ThreeVector vertex( 10.*myTrack.vfcalEntrance().vertex().X(),
223  10.*myTrack.vfcalEntrance().vertex().Y(),
224  10.*myTrack.vfcalEntrance().vertex().Z()+delZv); // in [mm]
225 
226  G4ThreeVector direction(myTrack.vfcalEntrance().Vect().X(),
227  myTrack.vfcalEntrance().Vect().Y(),
228  myTrack.vfcalEntrance().Vect().Z());
229 
230  bool ok;
231  double weight = 1.0; // rad. damage
232  int parCode = myTrack.type();
233  std::vector<FastHFShowerLibrary::Hit> hits =
234  getHits(vertex, direction, parCode, eGen, ok, weight, false);
235 
236  for (unsigned int i=0; i<hits.size(); ++i) {
237  G4ThreeVector pos = hits[i].position;
238  int depth = hits[i].depth;
239  double time = hits[i].time;
240 
241  if (isItinFidVolume (pos)) {
242  int det = 5;
243  int lay = 1;
244  uint32_t id = 0;
245  HcalNumberingFromDDD::HcalID tmp = numberingFromDDD->unitID(det, pos, depth, lay);
246  id = numberingScheme->getUnitID(tmp);
247 
248  CaloHitID current_id(id,time,myTrack.id());
249  std::map<CaloHitID,float>::iterator cellitr;
250  cellitr = hitMap.find(current_id);
251  if(cellitr==hitMap.end()) {
252  hitMap.insert(std::pair<CaloHitID,float>(current_id,1.0));
253  } else {
254  cellitr->second += 1.0;
255  }
256  } // end of isItinFidVolume check
257 
258  } // end loop over hits
259 
260 }
261 
262 bool FastHFShowerLibrary::isItinFidVolume (G4ThreeVector& hitPoint) {
263  bool flag = true;
264  if (applyFidCut) {
265  int npmt = HFFibreFiducial:: PMTNumber(hitPoint);
266 #ifdef DebugLog
267  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary::isItinFidVolume:#PMT= "
268  << npmt << " for hit point " << hitPoint;
269 #endif
270  if (npmt <= 0) flag = false;
271  }
272 #ifdef DebugLog
273  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary::isItinFidVolume: point "
274  << hitPoint << " return flag " << flag;
275 #endif
276  return flag;
277 }
278 
279 std::vector<FastHFShowerLibrary::Hit> FastHFShowerLibrary::getHits(const G4ThreeVector & hitPoint,
280  const G4ThreeVector & momDir, int parCode, double pin,
281  bool & ok, double weight, bool onlyLong) {
282 
283  std::vector<FastHFShowerLibrary::Hit> hit;
284  ok = false;
285  if (parCode == pi0PDG || parCode == etaPDG || parCode == nuePDG ||
286  parCode == numuPDG || parCode == nutauPDG || parCode == anuePDG ||
287  parCode == anumuPDG || parCode == anutauPDG || parCode == geantinoPDG)
288  return hit;
289  ok = true;
290 
291  double tSlice = 0.1*hitPoint.mag()/29.98;
292  double pz = momDir.z();
293  double zint = hitPoint.z();
294 
295 // if particle moves from interaction point or "backwards (halo)
296  int backward = 0;
297  if (pz * zint < 0.) backward = 1;
298 
299  double sphi = sin(momDir.phi());
300  double cphi = cos(momDir.phi());
301  double ctheta = cos(momDir.theta());
302  double stheta = sin(momDir.theta());
303 
304 #ifdef DebugLog
305  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: getHits " << parCode
306  << " of energy " << pin/GeV << " GeV"
307  << " dir.orts " << momDir.x() << ", " <<momDir.y()
308  << ", " << momDir.z() << " Pos x,y,z = "
309  << hitPoint.x() << "," << hitPoint.y() << ","
310  << hitPoint.z() << ","
311  << " sphi,cphi,stheta,ctheta = " << sphi
312  << "," << cphi << ", " << stheta << "," << ctheta;
313 #endif
314 
315  if (parCode == emPDG || parCode == epPDG || parCode == gammaPDG ) {
316  if (pin<pmom[nMomBin-1]) {
317  interpolate(0, pin);
318  } else {
319  extrapolate(0, pin);
320  }
321  } else {
322  if (pin<pmom[nMomBin-1]) {
323  interpolate(1, pin);
324  } else {
325  extrapolate(1, pin);
326  }
327  }
328 
329  int nHit = 0;
330  FastHFShowerLibrary::Hit oneHit;
331  for (int i = 0; i < npe; i++) {
332  double zv = std::abs(pe[i].z()); // abs local z
333 
334 #ifdef DebugLog
335  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: Hit " << i << " " << pe[i] << " zv " << zv;
336 #endif
337 
338  if (zv <= gpar[1] && pe[i].lambda() > 0 &&
339  (pe[i].z() >= 0 || (zv > gpar[0] && (!onlyLong)))) {
340  int depth = 1;
341  if (onlyLong) {
342  } else if (backward == 0) { // fully valid only for "front" particles
343  if (pe[i].z() < 0) depth = 2;// with "front"-simulated shower lib.
344  } else { // for "backward" particles - almost equal
345  double r = G4UniformRand(); // share between L and S fibers
346  if (r > 0.5) depth = 2;
347  }
348 
349 
350  // Updated coordinate transformation from local
351  // back to global using two Euler angles: phi and theta
352  double pex = pe[i].x();
353  double pey = pe[i].y();
354 
355  double xx = pex*ctheta*cphi - pey*sphi + zv*stheta*cphi;
356  double yy = pex*ctheta*sphi + pey*cphi + zv*stheta*sphi;
357  double zz = -pex*stheta + zv*ctheta;
358 
359  G4ThreeVector pos = hitPoint + G4ThreeVector(xx,yy,zz);
360  zv = std::abs(pos.z()) - gpar[4] - 0.5*gpar[1];
361  G4ThreeVector lpos = G4ThreeVector(pos.x(),pos.y(),zv);
362 
363  zv = fibre->zShift(lpos,depth,0); // distance to PMT !
364 
365  double r = pos.perp();
366  double p = fibre->attLength(pe[i].lambda());
367  double fi = pos.phi();
368  if (fi < 0) fi += twopi;
369  int isect = int(fi/dphi) + 1;
370  isect = (isect + 1) / 2;
371  double dfi = ((isect*2-1)*dphi - fi);
372  if (dfi < 0) dfi = -dfi;
373  double dfir = r * sin(dfi);
374 
375 #ifdef DebugLog
376  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: Position shift " << xx
377  << ", " << yy << ", " << zz << ": " << pos
378  << " R " << r << " Phi " << fi << " Section "
379  << isect << " R*Dfi " << dfir << " Dist " << zv;
380 #endif
381 
382  zz = std::abs(pos.z());
383  double r1 = G4UniformRand();
384  double r2 = G4UniformRand();
385  double r3 = -9999.;
386  if (backward) r3 = G4UniformRand();
387  if (!applyFidCut) dfir += gpar[5];
388 
389 #ifdef DebugLog
390  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: rLimits " << rInside(r)
391  << " attenuation " << r1 <<":" << exp(-p*zv)
392  << " r2 " << r2 << " r3 " << r3 << " rDfi "
393  << gpar[5] << " zz "
394  << zz << " zLim " << gpar[4] << ":"
395  << gpar[4]+gpar[1] << "\n"
396  << " rInside(r) :" << rInside(r)
397  << " r1 <= exp(-p*zv) :" << (r1 <= exp(-p*zv))
398  << " r2 <= probMax :" << (r2 <= probMax*weight)
399  << " r3 <= backProb :" << (r3 <= backProb)
400  << " dfir > gpar[5] :" << (dfir > gpar[5])
401  << " zz >= gpar[4] :" << (zz >= gpar[4])
402  << " zz <= gpar[4]+gpar[1] :"
403  << (zz <= gpar[4]+gpar[1]);
404 #endif
405 
406  if (rInside(r) && r1 <= exp(-p*zv) && r2 <= probMax*weight &&
407  dfir > gpar[5] && zz >= gpar[4] && zz <= gpar[4]+gpar[1] &&
408  r3 <= backProb && (depth != 2 || zz >= gpar[4]+gpar[0])) {
409 
410  oneHit.position = pos;
411  oneHit.depth = depth;
412  oneHit.time = (tSlice+(pe[i].t())+(fibre->tShift(lpos,depth,1)));
413  hit.push_back(oneHit);
414 #ifdef DebugLog
415  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: Final Hit " << nHit
416  <<" position " << (hit[nHit].position)
417  << " Depth " << (hit[nHit].depth) <<" Time "
418  << tSlice << ":" << pe[i].t() << ":"
419  << fibre->tShift(lpos,depth,1) << ":"
420  << (hit[nHit].time);
421 #endif
422  nHit++;
423  }
424 #ifdef DebugLog
425  else LogDebug("FastCalorimetry") << "HFShowerLibrary: REJECTED !!!";
426 #endif
427  if (onlyLong && zz >= gpar[4]+gpar[0] && zz <= gpar[4]+gpar[1]) {
428  r1 = G4UniformRand();
429  r2 = G4UniformRand();
430  if (rInside(r) && r1 <= exp(-p*zv) && r2 <= probMax && dfir > gpar[5]){
431  oneHit.position = pos;
432  oneHit.depth = 2;
433  oneHit.time = (tSlice+(pe[i].t())+(fibre->tShift(lpos,2,1)));
434  hit.push_back(oneHit);
435 #ifdef DebugLog
436  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: Final Hit " << nHit
437  << " position " << (hit[nHit].position)
438  << " Depth " << (hit[nHit].depth) <<" Time "
439  << tSlice << ":" << pe[i].t() << ":"
440  << fibre->tShift(lpos,2,1) << ":"
441  << (hit[nHit].time);
442 #endif
443  nHit++;
444  }
445  }
446  }
447  }
448 
449 #ifdef DebugLog
450  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary: Total Hits " << nHit
451  << " out of " << npe << " PE";
452 #endif
453 
454  if (nHit > npe && !onlyLong)
455  edm::LogWarning("FastCalorimetry") << "HFShowerLibrary: Hit buffer " << npe
456  << " smaller than " << nHit << " Hits";
457  return hit;
458 
459 }
460 
461 bool FastHFShowerLibrary::rInside(double r) {
462 
463  if (r >= rMin && r <= rMax) return true;
464  else return false;
465 }
466 
467 void FastHFShowerLibrary::getRecord(int type, int record) {
468 
469  int nrc = record-1;
470  photon.clear();
471  photo->clear();
472  if (type > 0) {
473  if (newForm) {
474  hadBranch->SetAddress(&photo);
475  hadBranch->GetEntry(nrc+totEvents);
476  } else {
477  hadBranch->SetAddress(&photon);
478  hadBranch->GetEntry(nrc);
479  }
480  } else {
481  if (newForm) {
482  emBranch->SetAddress(&photo);
483  } else {
484  emBranch->SetAddress(&photon);
485  }
486  emBranch->GetEntry(nrc);
487  }
488 
489 #ifdef DebugLog
490  int nPhoton = (newForm) ? photo->size() : photon.size();
491  LogDebug("FastCalorimetry") << "HFShowerLibrary::getRecord: Record " << record
492  << " of type " << type << " with " << nPhoton
493  << " photons";
494  for (int j = 0; j < nPhoton; j++)
495  if (newForm) LogDebug("FastCalorimetry") << "Photon " << j << " " << photo->at(j);
496  else LogDebug("FastCalorimetry") << "Photon " << j << " " << photon[j];
497 #endif
498 }
499 
500 void FastHFShowerLibrary::loadEventInfo(TBranch* branch) {
501 
502  if (branch) {
503  std::vector<HFShowerLibraryEventInfo> eventInfoCollection;
504  branch->SetAddress(&eventInfoCollection);
505  branch->GetEntry(0);
506  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary::loadEventInfo loads "
507  << " EventInfo Collection of size "
508  << eventInfoCollection.size() << " records";
509  totEvents = eventInfoCollection[0].totalEvents();
510  nMomBin = eventInfoCollection[0].numberOfBins();
511  evtPerBin = eventInfoCollection[0].eventsPerBin();
512  libVers = eventInfoCollection[0].showerLibraryVersion();
513  listVersion = eventInfoCollection[0].physListVersion();
514  pmom = eventInfoCollection[0].energyBins();
515  } else {
516  edm::LogInfo("FastCalorimetry") << "HFShowerLibrary::loadEventInfo loads "
517  << " EventInfo from hardwired numbers";
518  nMomBin = 16;
519  evtPerBin = 5000;
520  totEvents = nMomBin*evtPerBin;
521  libVers = 1.1;
522  listVersion = 3.6;
523  pmom = {2,3,5,7,10,15,20,30,50,75,100,150,250,350,500,1000};
524  }
525  for (int i=0; i<nMomBin; i++)
526  pmom[i] *= GeV;
527 }
528 
529 void FastHFShowerLibrary::interpolate(int type, double pin) {
530 
531 #ifdef DebugLog
532  LogDebug("FastCalorimetry") << "HFShowerLibrary:: Interpolate for Energy " <<pin/GeV
533  << " GeV with " << nMomBin << " momentum bins and "
534  << evtPerBin << " entries/bin -- total " << totEvents;
535 #endif
536 
537  int irc[2];
538  double w = 0.;
539  double r = G4UniformRand();
540 
541  if (pin<pmom[0]) {
542  w = pin/pmom[0];
543  irc[1] = int(evtPerBin*r) + 1;
544  irc[0] = 0;
545  } else {
546  for (int j=0; j<nMomBin-1; j++) {
547  if (pin >= pmom[j] && pin < pmom[j+1]) {
548  w = (pin-pmom[j])/(pmom[j+1]-pmom[j]);
549  if (j == nMomBin-2) {
550  irc[1] = int(evtPerBin*0.5*r);
551  } else {
552  irc[1] = int(evtPerBin*r);
553  }
554  irc[1] += (j+1)*evtPerBin + 1;
555  r = G4UniformRand();
556  irc[0] = int(evtPerBin*r) + 1 + j*evtPerBin;
557  if (irc[0]<0) {
558  edm::LogWarning("FastCalorimetry") << "HFShowerLibrary:: Illegal irc[0] = "
559  << irc[0] << " now set to 0";
560  irc[0] = 0;
561  } else if (irc[0] > totEvents) {
562  edm::LogWarning("FastCalorimetry") << "HFShowerLibrary:: Illegal irc[0] = "
563  << irc[0] << " now set to "<< totEvents;
564  irc[0] = totEvents;
565  }
566  }
567  }
568  }
569  if (irc[1]<1) {
570  edm::LogWarning("FastCalorimetry") << "HFShowerLibrary:: Illegal irc[1] = "
571  << irc[1] << " now set to 1";
572 
573  irc[1] = 1;
574  } else if (irc[1] > totEvents) {
575  edm::LogWarning("FastCalorimetry") << "HFShowerLibrary:: Illegal irc[1] = "
576  << irc[1] << " now set to "<< totEvents;
577 
578  irc[1] = totEvents;
579  }
580 
581 #ifdef DebugLog
582  LogDebug("FastCalorimetry") << "HFShowerLibrary:: Select records " << irc[0]
583  << " and " << irc[1] << " with weights " << 1-w
584  << " and " << w;
585 #endif
586 
587  pe.clear();
588  npe = 0;
589  int npold = 0;
590  for (int ir=0; ir < 2; ir++) {
591  if (irc[ir]>0) {
592  getRecord (type, irc[ir]);
593  int nPhoton = (newForm) ? photo->size() : photon.size();
594  npold += nPhoton;
595  for (int j=0; j<nPhoton; j++) {
596  r = G4UniformRand();
597  if ((ir==0 && r > w) || (ir > 0 && r < w)) {
598  storePhoton (j);
599  }
600  }
601  }
602  }
603 
604  if ((npe > npold || (npold == 0 && irc[0] > 0)) && !(npe == 0 && npold == 0))
605  edm::LogWarning("FastCalorimetry") << "HFShowerLibrary: Interpolation Warning =="
606  << " records " << irc[0] << " and " << irc[1]
607  << " gives a buffer of " << npold
608  << " photons and fills " << npe << " *****";
609 #ifdef DebugLog
610  else
611  LogDebug("FastCalorimetry") << "HFShowerLibrary: Interpolation == records "
612  << irc[0] << " and " << irc[1] << " gives a "
613  << "buffer of " << npold << " photons and fills "
614  << npe << " PE";
615  for (int j=0; j<npe; j++)
616  LogDebug("FastCalorimetry") << "Photon " << j << " " << pe[j];
617 #endif
618 
619 }
620 
621 void FastHFShowerLibrary::extrapolate(int type, double pin) {
622 
623  int nrec = int(pin/pmom[nMomBin-1]);
624  double w = (pin - pmom[nMomBin-1]*nrec)/pmom[nMomBin-1];
625  nrec++;
626 #ifdef DebugLog
627  LogDebug("FastCalorimetry") << "HFShowerLibrary:: Extrapolate for Energy " << pin
628  << " GeV with " << nMomBin << " momentum bins and "
629  << evtPerBin << " entries/bin -- total " << totEvents
630  << " using " << nrec << " records";
631 #endif
632 
633  std::vector<int> irc(nrec);
634 
635  for (int ir=0; ir<nrec; ir++) {
636  double r = G4UniformRand();
637  irc[ir] = int(evtPerBin*0.5*r) +(nMomBin-1)*evtPerBin + 1;
638  if (irc[ir]<1) {
639  edm::LogWarning("FastCalorimetry") << "HFShowerLibrary:: Illegal irc[" << ir
640  << "] = " << irc[ir] << " now set to 1";
641  irc[ir] = 1;
642  } else if (irc[ir] > totEvents) {
643  edm::LogWarning("FastCalorimetry") << "HFShowerLibrary:: Illegal irc[" << ir
644  << "] = " << irc[ir] << " now set to "
645  << totEvents;
646  irc[ir] = totEvents;
647 #ifdef DebugLog
648  } else {
649  LogDebug("FastCalorimetry") << "HFShowerLibrary::Extrapolation use irc["
650  << ir << "] = " << irc[ir];
651 #endif
652  }
653  }
654 
655  pe.clear();
656  npe = 0;
657  int npold = 0;
658  for (int ir=0; ir<nrec; ir++) {
659  if (irc[ir]>0) {
660  getRecord (type, irc[ir]);
661  int nPhoton = (newForm) ? photo->size() : photon.size();
662  npold += nPhoton;
663  for (int j=0; j<nPhoton; j++) {
664  double r = G4UniformRand();
665  if (ir != nrec-1 || r < w) {
666  storePhoton (j);
667  }
668  }
669 #ifdef DebugLog
670  LogDebug("FastCalorimetry") << "HFShowerLibrary: Record [" << ir << "] = "
671  << irc[ir] << " npold = " << npold;
672 #endif
673  }
674  }
675 #ifdef DebugLog
676  LogDebug("FastCalorimetry") << "HFShowerLibrary:: uses " << npold << " photons";
677 #endif
678 
679  if (npe > npold || npold == 0)
680  edm::LogWarning("FastCalorimetry") << "HFShowerLibrary: Extrapolation Warning == "
681  << nrec << " records " << irc[0] << ", "
682  << irc[1] << ", ... gives a buffer of " <<npold
683  << " photons and fills " << npe
684  << " *****";
685 #ifdef DebugLog
686  else
687  LogDebug("FastCalorimetry") << "HFShowerLibrary: Extrapolation == " << nrec
688  << " records " << irc[0] << ", " << irc[1]
689  << ", ... gives a buffer of " << npold
690  << " photons and fills " << npe << " PE";
691  for (int j=0; j<npe; j++)
692  LogDebug("FastCalorimetry") << "Photon " << j << " " << pe[j];
693 #endif
694 }
695 
696 void FastHFShowerLibrary::storePhoton(int j) {
697 
698  if (newForm) pe.push_back(photo->at(j));
699  else pe.push_back(photon[j]);
700 #ifdef DebugLog
701  LogDebug("FastCalorimetry") << "HFShowerLibrary: storePhoton " << j << " npe "
702  << npe << " " << pe[npe];
703 #endif
704 
705  npe++;
706 }
707 
708 std::vector<double> FastHFShowerLibrary::getDDDArray(const std::string & str,
709  const DDsvalues_type & sv,
710  int & nmin) {
711 
712 #ifdef DebugLog
713  LogDebug("FastCalorimetry") << "HFShowerLibrary:getDDDArray called for " << str
714  << " with nMin " << nmin;
715 #endif
716 
717  DDValue value(str);
718  if (DDfetch(&sv,value)) {
719 #ifdef DebugLog
720  LogDebug("FastCalorimetry") << "HFShowerLibrary:getDDDArray value " << value;
721 #endif
722 
723  const std::vector<double> & fvec = value.doubles();
724  int nval = fvec.size();
725  if (nmin > 0) {
726  if (nval < nmin) {
727  edm::LogError("FastCalorimetry") << "HFShowerLibrary : # of " << str
728  << " bins " << nval << " < " << nmin
729  << " ==> illegal";
730  throw cms::Exception("Unknown", "HFShowerLibrary")
731  << "nval < nmin for array " << str << "\n";
732  }
733  } else {
734  if (nval < 2) {
735  edm::LogError("FastCalorimetry") << "HFShowerLibrary : # of " << str
736  << " bins " << nval << " < 2 ==> illegal"
737  << " (nmin=" << nmin << ")";
738  throw cms::Exception("Unknown", "HFShowerLibrary")
739  << "nval < 2 for array " << str << "\n";
740  }
741  }
742  nmin = nval;
743 
744  return fvec;
745  } else {
746  edm::LogError("FastCalorimetry") << "HFShowerLibrary : cannot get array " << str;
747 
748  throw cms::Exception("Unknown", "HFShowerLibrary")
749  << "cannot get array " << str << "\n";
750  }
751 }
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:501
type
type
Definition: HCALResponse.h:21
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
FastHFShowerLibrary::numberingScheme
HcalNumberingScheme * numberingScheme
Definition: FastHFShowerLibrary.h:81
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T 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
FastHFShowerLibrary::getDDDArray
std::vector< double > getDDDArray(const std::string &, const DDsvalues_type &, int &)
Definition: FastHFShowerLibrary.cc:708
FSimTrack::vfcalEntrance
const RawParticle & vfcalEntrance() const
The particle at VFCAL entrance.
Definition: FSimTrack.h:139
w
const double w
Definition: UKUtility.cc:23
record
JetCorrectorParameters::Record record
Definition: classes.h:7
FastHFShowerLibrary::recoHFShowerLibrary
void recoHFShowerLibrary(const FSimTrack &myTrack)
Definition: FastHFShowerLibrary.cc:207
HcalNumberingScheme::getUnitID
virtual uint32_t getUnitID(const HcalNumberingFromDDD::HcalID &id)
Definition: HcalNumberingScheme.cc:22
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
FastHFShowerLibrary::pe
HFShowerPhotonCollection pe
Definition: FastHFShowerLibrary.h:100
alcazmumu_cfi.filter
tuple filter
Definition: alcazmumu_cfi.py:34
FastHFShowerLibrary::loadEventInfo
void loadEventInfo(TBranch *)
Definition: FastHFShowerLibrary.cc:500
DDfetch
bool DDfetch(const DDsvalues_type *, DDValue &)
helper for retrieving DDValues from DDsvalues_type *.
Definition: DDsvalues.cc:80
detailsBasic3DVector::z
float float float z
Definition: extBasic3DVector.h:15
HFFibreFiducial::PMTNumber
static int PMTNumber(const G4ThreeVector &pe_effect)
Definition: HFFibreFiducial.cc:12
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
numberingScheme
static TrackerG4SimHitNumberingScheme & numberingScheme(const DDCompactView &cpv, const GeometricDet &det)
Definition: TkAccumulatingSensitiveDetector.cc:57
FastHFShowerLibrary::numberingFromDDD
HcalNumberingFromDDD * numberingFromDDD
Definition: FastHFShowerLibrary.h:82
FastHFShowerLibrary::interpolate
void interpolate(int, double)
Definition: FastHFShowerLibrary.cc:529
HFShowerPhotonCollection
std::vector< HFShowerPhoton > HFShowerPhotonCollection
Definition: HFShowerPhoton.h:40
FastHFFibre::tShift
double tShift(const G4ThreeVector &point, int depth, int fromEndAbs=0)
Definition: FastHFFibre.cc:128
HLT_25ns14e33_v1_cff.nEta
tuple nEta
Definition: HLT_25ns14e33_v1_cff.py:1119
FastHFShowerLibrary::getHits
std::vector< Hit > getHits(const G4ThreeVector &p, const G4ThreeVector &v, int parCode, double parEnergy, bool &ok, double weight, bool onlyLong=false)
Definition: FastHFShowerLibrary.cc:279
FastHFShowerLibrary::extrapolate
void extrapolate(int, double)
Definition: FastHFShowerLibrary.cc:621
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
RawParticle::Y
double Y() const
y of vertex
Definition: RawParticle.h:275
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
FastHFFibre::attLength
double attLength(double lambda)
Definition: FastHFFibre.cc:110
RawParticle::Z
double Z() const
z of vertex
Definition: RawParticle.h:276
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
j
int j
Definition: DBlmapReader.cc:9
FSimTrack::onVFcal
int onVFcal() const
Definition: FSimTrack.h:111
FastHFShowerLibrary::photo
HFShowerPhotonCollection * photo
Definition: FastHFShowerLibrary.h:101
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
FastHFShowerLibrary::pmom
std::vector< double > pmom
Definition: FastHFShowerLibrary.h:89
HLT_25ns14e33_v1_cff.depth
tuple depth
Definition: HLT_25ns14e33_v1_cff.py:6940
RawParticle::vertex
const XYZTLorentzVector & vertex() const
the vertex fourvector
Definition: RawParticle.h:285
FastHFShowerLibrary::photon
HFShowerPhotonCollection photon
Definition: FastHFShowerLibrary.h:102
FastHFShowerLibrary::FastHFShowerLibrary
FastHFShowerLibrary(edm::ParameterSet const &p)
Definition: FastHFShowerLibrary.cc:28
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
DDFilteredView::mergedSpecifics
DDsvalues_type mergedSpecifics() const
Definition: DDFilteredView.cc:66
RawParticle::X
double X() const
x of vertex
Definition: RawParticle.h:274
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
FastHFShowerLibrary::gpar
std::vector< double > gpar
Definition: FastHFShowerLibrary.h:93
tmp
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
HLT_25ns14e33_v1_cff.nPhi
tuple nPhi
Definition: HLT_25ns14e33_v1_cff.py:1120
CoreSimTrack::type
int type() const
particle type (HEP PDT convension)
Definition: CoreSimTrack.h:25
FastHFFibre::zShift
double zShift(const G4ThreeVector &point, int depth, int fromEndAbs=0)
Definition: FastHFFibre.cc:140
FSimTrack::id
int id() const
the index in FBaseSimEvent and other vectors
Definition: FSimTrack.h:86
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
FastHFShowerLibrary::initHFShowerLibrary
void const initHFShowerLibrary(const edm::EventSetup &)
Definition: FastHFShowerLibrary.cc:131
HcalNumberingFromDDD::unitID
HcalID unitID(int det, const CLHEP::Hep3Vector &pos, int depth, int lay=-1) const
Definition: HcalNumberingFromDDD.cc:31
FastHFShowerLibrary::hitMap
std::map< CaloHitID, float > hitMap
Definition: FastHFShowerLibrary.h:83
FastHFShowerLibrary::isItinFidVolume
bool isItinFidVolume(G4ThreeVector &)
Definition: FastHFShowerLibrary.cc:262
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