CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CalorimetryManager.cc
Go to the documentation of this file.
1 //updated by Reza Goldouzian
2 //Framework headers
5 
6 // Fast Simulation headers
18 //#include "FastSimulation/Utilities/interface/Histos.h"
28 // New headers for Muon Mip Simulation
31 // Muon Brem
33 
34 //Gflash Hadronic Model
42 
43 //FastHFShowerLibrary
45 
46 // STL headers
47 #include <vector>
48 #include <iostream>
49 
50 //CMSSW headers
53 //#include "DataFormats/EcalDetId/interface/EcalDetId.h"
55 
58 
59 //ROOT headers
60 #include "TROOT.h"
61 #include "TH1.h"
62 
63 using namespace edm;
64 
67 
68 std::vector<std::pair<int, float> > CalorimetryManager::myZero_ =
69  std::vector<std::pair<int, float> >(1, std::pair<int, float>(0, 0.));
70 
71 CalorimetryManager::CalorimetryManager() : myCalorimeter_(nullptr), initialized_(false) { ; }
72 
74  const edm::ParameterSet& fastCalo,
75  const edm::ParameterSet& fastMuECAL,
76  const edm::ParameterSet& fastMuHCAL,
77  const edm::ParameterSet& parGflash)
78  : mySimEvent(aSimEvent),
85 
86  //Gflash
87  theProfile = new GflashHadronShowerProfile(parGflash);
88  thePiKProfile = new GflashPiKShowerProfile(parGflash);
91 
92  // FastHFShowerLibrary
94 
95  readParameters(fastCalo);
96 
97  myCalorimeter_ = new CaloGeometryHelper(fastCalo);
98  myHDResponse_ = new HCALResponse(fastCalo.getParameter<edm::ParameterSet>("HCALResponse"));
99  myHSParameters_ = new HSParameters(fastCalo.getParameter<edm::ParameterSet>("HSParameters"));
100 
101  // Material Effects for Muons in ECAL (only EnergyLoss implemented so far)
102  if (fastMuECAL.getParameter<bool>("PairProduction") || fastMuECAL.getParameter<bool>("Bremsstrahlung") ||
103  fastMuECAL.getParameter<bool>("MuonBremsstrahlung") || fastMuECAL.getParameter<bool>("EnergyLoss") ||
104  fastMuECAL.getParameter<bool>("MultipleScattering"))
105  theMuonEcalEffects = new MaterialEffects(fastMuECAL);
106 
107  // Material Effects for Muons in HCAL (only EnergyLoss implemented so far)
108  if (fastMuHCAL.getParameter<bool>("PairProduction") || fastMuHCAL.getParameter<bool>("Bremsstrahlung") ||
109  fastMuHCAL.getParameter<bool>("MuonBremsstrahlung") || fastMuHCAL.getParameter<bool>("EnergyLoss") ||
110  fastMuHCAL.getParameter<bool>("MultipleScattering"))
111  theMuonHcalEffects = new MaterialEffects(fastMuHCAL);
112 
113  if (fastCalo.exists("ECALResponseScaling")) {
114  ecalCorrection = std::unique_ptr<KKCorrectionFactors>(
115  new KKCorrectionFactors(fastCalo.getParameter<edm::ParameterSet>("ECALResponseScaling")));
116  }
117 }
118 
120  EBMapping_.clear();
121  EEMapping_.clear();
122  HMapping_.clear();
123  ESMapping_.clear();
124  muonSimTracks.clear();
125  savedMuonSimTracks.clear();
126 }
127 
129  if (myCalorimeter_)
130  delete myCalorimeter_;
131  if (myHDResponse_)
132  delete myHDResponse_;
133 
134  if (theMuonEcalEffects)
135  delete theMuonEcalEffects;
136  if (theMuonHcalEffects)
137  delete theMuonHcalEffects;
138 
139  if (theProfile)
140  delete theProfile;
141 
142  if (theHFShowerLibrary)
143  delete theHFShowerLibrary;
144 }
145 
147  if (!evtsToDebug_.empty()) {
148  std::vector<unsigned int>::const_iterator itcheck =
149  find(evtsToDebug_.begin(), evtsToDebug_.end(), mySimEvent->id().event());
150  debug_ = (itcheck != evtsToDebug_.end());
151  if (debug_)
152  mySimEvent->print();
153  }
154 
155  initialize(random);
156 
157  LogInfo("FastCalorimetry") << "Reconstructing " << (int)mySimEvent->nTracks() << " tracks." << std::endl;
158  for (int fsimi = 0; fsimi < (int)mySimEvent->nTracks(); ++fsimi) {
159  FSimTrack& myTrack = mySimEvent->track(fsimi);
160 
161  reconstructTrack(myTrack, random);
162  } // particle loop
163 
164 } // reconstruct
165 
167  // Clear the content of the calorimeters
168  if (!initialized_) {
169  theHFShowerLibrary->SetRandom(random);
170 
171  // Check if the preshower is really available
173  edm::LogWarning("CalorimetryManager")
174  << " WARNING: The preshower simulation has been turned on; but no preshower geometry is available "
175  << std::endl;
176  edm::LogWarning("CalorimetryManager") << " Disabling the preshower simulation " << std::endl;
177  simulatePreshower_ = false;
178  }
179 
180  initialized_ = true;
181  }
182  clean();
183 }
184 
186  int pid = abs(myTrack.type());
187 
188  if (debug_) {
189  LogInfo("FastCalorimetry") << " ===> pid = " << pid << std::endl;
190  }
191 
192  // Check that the particle hasn't decayed
193  if (myTrack.noEndVertex()) {
194  // Simulate energy smearing for photon and electrons
195  float charge_ = (float)(myTrack.charge());
196  if (pid == 11 || pid == 22) {
197  if (myTrack.onEcal())
198  EMShowerSimulation(myTrack, random);
199  else if (myTrack.onVFcal()) {
200  if (useShowerLibrary) {
202  myHDResponse_->correctHF(myTrack.hcalEntrance().e(), abs(myTrack.type()));
204  } else
205  reconstructHCAL(myTrack, random);
206  }
207  } // electron or photon
208  else if (pid == 13 || pid == 1000024 || (pid > 1000100 && pid < 1999999 && fabs(charge_) > 0.001)) {
209  MuonMipSimulation(myTrack, random);
210  }
211  // Simulate energy smearing for hadrons (i.e., everything
212  // but muons... and SUSY particles that deserve a special
213  // treatment.
214  else if (pid < 1000000) {
215  if (myTrack.onHcal() || myTrack.onVFcal()) {
216  if (optionHDSim_ == 0)
217  reconstructHCAL(myTrack, random);
218  else
219  HDShowerSimulation(myTrack, random);
220  }
221  } // pid < 1000000
222  } // myTrack.noEndVertex()
223 }
224 
225 // Simulation of electromagnetic showers in PS, ECAL, HCAL
227  std::vector<const RawParticle*> thePart;
228  double X0depth;
229 
230  if (debug_) {
231  LogInfo("FastCalorimetry") << " EMShowerSimulation " << myTrack << std::endl;
232  }
233 
234  // The Particle at ECAL entrance
235  myPart = myTrack.ecalEntrance();
236 
237  // protection against infinite loop.
238  if (myTrack.type() == 22 && myPart.e() < 0.055)
239  return;
240 
241  // Barrel or Endcap ?
242  int onEcal = myTrack.onEcal();
243  int onHcal = myTrack.onHcal();
244  int onLayer1 = myTrack.onLayer1();
245  int onLayer2 = myTrack.onLayer2();
246 
247  // The entrance in ECAL
248  XYZPoint ecalentrance = myPart.vertex().Vect();
249 
250  // The preshower
251  PreshowerHitMaker* myPreshower = nullptr;
252  if (simulatePreshower_ && (onLayer1 || onLayer2)) {
253  XYZPoint layer1entrance, layer2entrance;
255  if (onLayer1) {
256  layer1entrance = XYZPoint(myTrack.layer1Entrance().vertex().Vect());
257  dir1 = XYZVector(myTrack.layer1Entrance().Vect().Unit());
258  }
259  if (onLayer2) {
260  layer2entrance = XYZPoint(myTrack.layer2Entrance().vertex().Vect());
261  dir2 = XYZVector(myTrack.layer2Entrance().Vect().Unit());
262  }
263  myPreshower =
264  new PreshowerHitMaker(myCalorimeter_, layer1entrance, dir1, layer2entrance, dir2, aLandauGenerator, random);
265  myPreshower->setMipEnergy(mipValues_[0], mipValues_[1]);
266  }
267 
268  // The ECAL Properties
271  myCalorimeter_->layer1Properties(onLayer1),
272  myCalorimeter_->layer2Properties(onLayer2),
275  RCFactor_,
276  RTFactor_);
277 
278  // Photons : create an e+e- pair
279  if (myTrack.type() == 22) {
280  // Depth for the first e+e- pair creation (in X0)
281  X0depth = -log(random->flatShoot()) * (9. / 7.);
282 
283  // Initialization
284  double eMass = 0.000510998902;
285  double xe = 0;
286  double xm = eMass / myPart.e();
287  double weight = 0.;
288 
289  // Generate electron energy between emass and eGamma-emass
290  do {
291  xe = random->flatShoot() * (1. - 2. * xm) + xm;
292  weight = 1. - 4. / 3. * xe * (1. - xe);
293  } while (weight < random->flatShoot());
294 
295  // Protection agains infinite loop in Famos Shower
296  if (myPart.e() * xe < 0.055 || myPart.e() * (1. - xe) < 0.055) {
297  if (myPart.e() > 0.055)
298  thePart.push_back(&myPart);
299 
300  } else {
301  myElec = (myPart.momentum()) * xe;
302  myPosi = (myPart.momentum()) * (1. - xe);
305  thePart.push_back(&myElec);
306  thePart.push_back(&myPosi);
307  }
308  // Electrons
309  } else {
310  X0depth = 0.;
311  if (myPart.e() > 0.055)
312  thePart.push_back(&myPart);
313  }
314 
315  // After the different protections, this shouldn't happen.
316  if (thePart.empty()) {
317  if (myPreshower == nullptr)
318  return;
319  delete myPreshower;
320  return;
321  }
322 
323  // find the most energetic particle
324  double maxEnergy = -1.;
325  for (unsigned ip = 0; ip < thePart.size(); ++ip)
326  if (thePart[ip]->e() > maxEnergy)
327  maxEnergy = thePart[ip]->e();
328 
329  // Initialize the Grid in ECAL
330  int size = gridSize_;
331  if (maxEnergy > 100)
332  size = 11;
333 
334  EMShower theShower(random, aGammaGenerator, &showerparam, &thePart, nullptr, nullptr, bFixedLength_);
335 
336  double maxShower = theShower.getMaximumOfShower();
337  if (maxShower > 20.)
338  maxShower = 2.; // simple pivot-searching protection
339 
340  double depth((X0depth + maxShower) * myCalorimeter_->ecalProperties(onEcal)->radLenIncm());
341  XYZPoint meanShower = ecalentrance + myPart.Vect().Unit() * depth;
342 
343  // The closest crystal
344  DetId pivot(myCalorimeter_->getClosestCell(meanShower, true, onEcal == 1));
345 
346  if (pivot.subdetId() == 0) { // further protection against avbsence of pivot
347  edm::LogWarning("CalorimetryManager")
348  << "Pivot for egamma e = " << myTrack.hcalEntrance().e() << " is not found at depth " << depth
349  << " and meanShower coordinates = " << meanShower << std::endl;
350  if (myPreshower)
351  delete myPreshower;
352  return;
353  }
354 
355  EcalHitMaker myGrid(myCalorimeter_, ecalentrance, pivot, onEcal, size, 0, random);
356  // ^^^^
357  // for EM showers
360 
361  //maximumdepth dependence of the radiusfactorbehindpreshower
362  //First tuning: Shilpi Jain (Mar-Apr 2010); changed after tuning - Feb-July - Shilpi Jain
363  /* **************
364  myGrid.setRadiusFactor(radiusFactor_);
365  if(onLayer1 || onLayer2)
366  {
367  float b = radiusPreshowerCorrections_[0];
368  float a = radiusFactor_*( 1.+radiusPreshowerCorrections_[1]*radiusPreshowerCorrections_[0] );
369  float maxdepth = X0depth+theShower.getMaximumOfShower();
370  float newRadiusFactor = radiusFactor_;
371  if(myPart.e()<=250.)
372  {
373  newRadiusFactor = a/(1.+b*maxdepth);
374  }
375  myGrid.setRadiusFactor(newRadiusFactor);
376  }
377  else // otherwise use the normal radius factor
378  {
379  myGrid.setRadiusFactor(radiusFactor_);
380  }
381  ************** */
382  if (myTrack.onEcal() == 2) // if on EE
383  {
384  if ((onLayer1 || onLayer2) && myPart.e() <= 250.) {
385  double maxdepth = X0depth + theShower.getMaximumOfShower();
386  double newRadiusFactor = radiusFactorEE_ * aTerm / (1. + bTerm * maxdepth);
387  myGrid.setRadiusFactor(newRadiusFactor);
388  } else // otherwise use the normal radius factor
389  {
391  }
392  } //if(myTrack.onEcal() == 2)
393  else // else if on EB
394  {
396  }
397  //(end of) changed after tuning - Feb-July - Shilpi Jain
398 
400 
401  // The shower simulation
402  myGrid.setTrackParameters(myPart.Vect().Unit(), X0depth, myTrack);
403 
404  if (myPreshower)
405  theShower.setPreshower(myPreshower);
406 
407  HcalHitMaker myHcalHitMaker(myGrid, (unsigned)0);
408 
409  theShower.setGrid(&myGrid);
410  theShower.setHcal(&myHcalHitMaker);
411  theShower.compute();
412 
413  // calculate the total simulated energy for this particle
414  float simE = 0;
415  for (const auto& mapIterator : myGrid.getHits()) {
416  simE += mapIterator.second;
417  }
418 
419  auto scale = ecalCorrection
420  ? ecalCorrection->getScale(myTrack.ecalEntrance().e(), std::abs(myTrack.ecalEntrance().eta()), simE)
421  : 1.;
422 
423  // Save the hits !
424  updateECAL(myGrid.getHits(), onEcal, myTrack.id(), scale);
425 
426  // Now fill the HCAL hits
427  updateHCAL(myHcalHitMaker.getHits(), myTrack.id());
428 
429  // delete the preshower
430  if (myPreshower != nullptr) {
431  updatePreshower(myPreshower->getHits(), myTrack.id());
432  delete myPreshower;
433  }
434 }
435 
437  int hit;
438  int pid = abs(myTrack.type());
439  if (debug_) {
440  LogInfo("FastCalorimetry") << " reconstructHCAL " << myTrack << std::endl;
441  }
442 
443  XYZTLorentzVector trackPosition;
444  if (myTrack.onHcal()) {
445  trackPosition = myTrack.hcalEntrance().vertex();
446  hit = myTrack.onHcal() - 1;
447  } else {
448  trackPosition = myTrack.vfcalEntrance().vertex();
449  hit = 2;
450  }
451 
452  double pathEta = trackPosition.eta();
453  double pathPhi = trackPosition.phi();
454 
455  double EGen = myTrack.hcalEntrance().e();
456  double emeas = 0.;
457 
458  float charge_ = (float)myTrack.charge();
459  if (pid == 13 || pid == 1000024 || (pid > 1000100 && pid < 1999999 && fabs(charge_) > 0.001)) {
460  emeas = myHDResponse_->responseHCAL(0, EGen, pathEta, 2, random); // 2=muon
461  if (debug_)
462  LogInfo("FastCalorimetry") << "CalorimetryManager::reconstructHCAL - MUON !!!" << std::endl;
463  } else if (pid == 22 || pid == 11) {
464  emeas = myHDResponse_->responseHCAL(0, EGen, pathEta, 0, random); // last par. = 0 = e/gamma
465  if (debug_)
466  LogInfo("FastCalorimetry") << "CalorimetryManager::reconstructHCAL - e/gamma !!!" << std::endl;
467  } else {
468  emeas = myHDResponse_->getHCALEnergyResponse(EGen, hit, random);
469  }
470 
471  if (debug_)
472  LogInfo("FastCalorimetry") << "CalorimetryManager::reconstructHCAL - on-calo "
473  << " eta = " << pathEta << " phi = " << pathPhi << " Egen = " << EGen
474  << " Emeas = " << emeas << std::endl;
475 
476  if (emeas > 0.) {
477  DetId cell = myCalorimeter_->getClosestCell(trackPosition.Vect(), false, false);
478  double tof =
479  (((HcalGeometry*)(myCalorimeter_->getHcalGeometry()))->getPosition(cell).mag()) / 29.98; //speed of light
480  CaloHitID current_id(cell.rawId(), tof, myTrack.id());
481  std::map<CaloHitID, float> hitMap;
482  hitMap[current_id] = emeas;
483  updateHCAL(hitMap, myTrack.id());
484  }
485 }
486 
488  // const edm::ParameterSet& fastCalo){
489 
490  theHFShowerLibrary->SetRandom(random);
491 
492  // TimeMe t(" FASTEnergyReconstructor::HDShower");
493  const XYZTLorentzVector& moment = myTrack.momentum();
494 
495  if (debug_)
496  LogInfo("FastCalorimetry") << "CalorimetryManager::HDShowerSimulation - track param." << std::endl
497  << " eta = " << moment.eta() << std::endl
498  << " phi = " << moment.phi() << std::endl
499  << " et = " << moment.Et() << std::endl
500  << " e = " << myTrack.hcalEntrance().e() << std::endl;
501 
502  if (debug_) {
503  LogInfo("FastCalorimetry") << " HDShowerSimulation " << myTrack << std::endl;
504  }
505 
506  int hit;
507 
508  XYZTLorentzVector trackPosition;
509  if (myTrack.onEcal()) {
510  trackPosition = myTrack.ecalEntrance().vertex();
511  hit = myTrack.onEcal() - 1; //
512  myPart = myTrack.ecalEntrance();
513  } else if (myTrack.onVFcal()) {
514  trackPosition = myTrack.vfcalEntrance().vertex();
515  hit = 2;
516  myPart = myTrack.vfcalEntrance();
517  } else {
518  LogInfo("FastCalorimetry") << " The particle is not in the acceptance " << std::endl;
519  return;
520  }
521 
522  // int onHCAL = hit + 1; - specially for myCalorimeter->hcalProperties(onHCAL)
523  // (below) to get VFcal properties ...
524  int onHCAL = hit + 1;
525  int onECAL = myTrack.onEcal();
526 
527  double pathEta = trackPosition.eta();
528  double pathPhi = trackPosition.phi();
529 
530  double eint = moment.e();
531  double eGen = myTrack.hcalEntrance().e();
532 
533  double emeas = 0.;
534  double pmip = myHDResponse_->getMIPfraction(eGen, pathEta);
535 
536  //===========================================================================
537  if (eGen > 0.) {
538  // ECAL and HCAL properties to get
539  HDShowerParametrization theHDShowerparam(
541 
542  //Making ECAL Grid (and segments calculation)
543  XYZPoint caloentrance;
544  XYZVector direction;
545  if (myTrack.onEcal()) {
546  caloentrance = myTrack.ecalEntrance().vertex().Vect();
547  direction = myTrack.ecalEntrance().Vect().Unit();
548  } else if (myTrack.onHcal()) {
549  caloentrance = myTrack.hcalEntrance().vertex().Vect();
550  direction = myTrack.hcalEntrance().Vect().Unit();
551  } else {
552  caloentrance = myTrack.vfcalEntrance().vertex().Vect();
553  direction = myTrack.vfcalEntrance().Vect().Unit();
554  }
555 
556  if (debug_)
557  LogInfo("FastCalorimetry") << "CalorimetryManager::HDShowerSimulation - on-calo 1 " << std::endl
558  << " onEcal = " << myTrack.onEcal() << std::endl
559  << " onHcal = " << myTrack.onHcal() << std::endl
560  << " onVFcal = " << myTrack.onVFcal() << std::endl
561  << " position = " << caloentrance << std::endl;
562 
563  DetId pivot;
564  if (myTrack.onEcal()) {
565  pivot = myCalorimeter_->getClosestCell(caloentrance, true, myTrack.onEcal() == 1);
566  } else if (myTrack.onHcal()) {
567  pivot = myCalorimeter_->getClosestCell(caloentrance, false, false);
568  }
569 
570  EcalHitMaker myGrid(
571  myCalorimeter_, caloentrance, pivot, pivot.null() ? 0 : myTrack.onEcal(), hdGridSize_, 1, random);
572  // 1=HAD shower
573 
574  myGrid.setTrackParameters(direction, 0, myTrack);
575  // Build the FAMOS HCAL
576  HcalHitMaker myHcalHitMaker(myGrid, (unsigned)1);
577 
578  // Shower simulation
579  bool status = false;
580  int mip = 2;
581  // Use HFShower for HF
582  if (!myTrack.onEcal() && !myTrack.onHcal()) {
583  // Warning : We give here the particle energy with the response
584  // but without the resolution/gaussian smearing
585  // For HF, the resolution is due to the PE statistic
586 
587  if (useShowerLibrary) {
589  status = true;
590  } else {
591  HFShower theShower(random, &theHDShowerparam, &myGrid, &myHcalHitMaker, onECAL, eGen);
592  // eGen);
593  // e); // PV Warning : temporarly set the energy to the generated E
594 
595  status = theShower.compute();
596  }
597  } else {
598  if (hdSimMethod_ == 0) {
599  HDShower theShower(random, &theHDShowerparam, &myGrid, &myHcalHitMaker, onECAL, eGen, pmip);
600  status = theShower.compute();
601  mip = theShower.getmip();
602  } else if (hdSimMethod_ == 1) {
603  HDRShower theShower(random, &theHDShowerparam, &myGrid, &myHcalHitMaker, onECAL, eGen);
604  status = theShower.computeShower();
605  mip = 2;
606  } else if (hdSimMethod_ == 2) {
607  //dynamically loading a corresponding profile by the particle type
608  int particleType = myTrack.type();
610  if (particleType == -2212)
612  else if (particleType == 2212)
614 
615  //input variables for GflashHadronShowerProfile
616  int showerType = 99 + myTrack.onEcal();
617  double globalTime = 150.0; // a temporary reference hit time in nanosecond
618  float charge = (float)(myTrack.charge());
619  Gflash3Vector gfpos(trackPosition.X(), trackPosition.Y(), trackPosition.Z());
620  Gflash3Vector gfmom(moment.X(), moment.Y(), moment.Z());
621 
622  theProfile->initialize(showerType, eGen, globalTime, charge, gfpos, gfmom);
625 
626  //make hits
627  std::vector<GflashHit>& gflashHitList = theProfile->getGflashHitList();
628  std::vector<GflashHit>::const_iterator spotIter = gflashHitList.begin();
629  std::vector<GflashHit>::const_iterator spotIterEnd = gflashHitList.end();
630 
632 
633  for (; spotIter != spotIterEnd; spotIter++) {
634  double pathLength = theProfile->getGflashShowino()->getPathLengthAtShower() +
635  (30 * 100 / eGen) * (spotIter->getTime() - globalTime);
636 
637  double currentDepth = std::max(0.0, pathLength - theProfile->getGflashShowino()->getPathLengthOnEcal());
638 
639  //find the the showino position at the currentDepth
640  GflashTrajectoryPoint trajectoryPoint;
641  theProfile->getGflashShowino()->getHelix()->getGflashTrajectoryPoint(trajectoryPoint, pathLength);
642  Gflash3Vector positionAtCurrentDepth = trajectoryPoint.getPosition();
643  //find radial distrance
644  Gflash3Vector lateralDisplacement = positionAtCurrentDepth - spotIter->getPosition() / CLHEP::cm;
645  double rShower = lateralDisplacement.r();
646  double azimuthalAngle = lateralDisplacement.phi();
647 
648  whichCalor = Gflash::getCalorimeterNumber(positionAtCurrentDepth);
649 
650  if (whichCalor == Gflash::kESPM || whichCalor == Gflash::kENCA) {
651  bool statusPad = myGrid.getPads(currentDepth, true);
652  if (!statusPad)
653  continue;
654  myGrid.setSpotEnergy(1.2 * spotIter->getEnergy() / CLHEP::GeV);
655  myGrid.addHit(rShower / Gflash::intLength[Gflash::kESPM], azimuthalAngle, 0);
656  } else if (whichCalor == Gflash::kHB || whichCalor == Gflash::kHE) {
657  bool setHDdepth = myHcalHitMaker.setDepth(currentDepth, true);
658  if (!setHDdepth)
659  continue;
660  myHcalHitMaker.setSpotEnergy(1.4 * spotIter->getEnergy() / CLHEP::GeV);
661  myHcalHitMaker.addHit(rShower / Gflash::intLength[Gflash::kHB], azimuthalAngle, 0);
662  }
663  }
664  status = true;
665  } else {
666  edm::LogInfo("FastSimulationCalorimetry") << " SimMethod " << hdSimMethod_ << " is NOT available ";
667  }
668  }
669 
670  if (status) {
671  // Here to switch between simple formulae and parameterized response
672  if (optionHDSim_ == 1) {
673  emeas = myHDResponse_->getHCALEnergyResponse(eGen, hit, random);
674  } else { // optionHDsim == 2
675  emeas = myHDResponse_->responseHCAL(mip, eGen, pathEta, 1, random); // 1=hadron
676  }
677 
678  double correction = emeas / eGen;
679 
680  // RespCorrP factors (ECAL and HCAL separately) calculation
681  respCorr(eint);
682 
683  if (debug_)
684  LogInfo("FastCalorimetry") << "CalorimetryManager::HDShowerSimulation - on-calo 2" << std::endl
685  << " eta = " << pathEta << std::endl
686  << " phi = " << pathPhi << std::endl
687  << " Egen = " << eGen << std::endl
688  << " Emeas = " << emeas << std::endl
689  << " corr = " << correction << std::endl
690  << " mip = " << mip << std::endl;
691 
692  if (myTrack.onEcal() > 0) {
693  // Save ECAL hits
694  updateECAL(myGrid.getHits(), onECAL, myTrack.id(), correction * ecorr);
695  }
696 
697  // Save HCAL hits
698  if (myTrack.onVFcal() && useShowerLibrary) {
699  myHDResponse_->correctHF(eGen, abs(myTrack.type()));
701  } else
702  updateHCAL(myHcalHitMaker.getHits(), myTrack.id(), correction * hcorr);
703 
704  } else { // shower simulation failed
705  if (myTrack.onHcal() || myTrack.onVFcal()) {
706  DetId cell = myCalorimeter_->getClosestCell(trackPosition.Vect(), false, false);
707  double tof =
708  (((HcalGeometry*)(myCalorimeter_->getHcalGeometry()))->getPosition(cell).mag()) / 29.98; //speed of light
709  CaloHitID current_id(cell.rawId(), tof, myTrack.id());
710  std::map<CaloHitID, float> hitMap;
711  hitMap[current_id] = emeas;
712  updateHCAL(hitMap, myTrack.id());
713  if (debug_)
714  LogInfo("FastCalorimetry") << " HCAL simple cell " << cell.rawId() << " added E = " << emeas << std::endl;
715  }
716  }
717 
718  } // e > 0. ...
719 
720  if (debug_)
721  LogInfo("FastCalorimetry") << std::endl << " FASTEnergyReconstructor::HDShowerSimulation finished " << std::endl;
722 }
723 
725  // TimeMe t(" FASTEnergyReconstructor::HDShower");
726  XYZTLorentzVector moment = myTrack.momentum();
727 
728  // Backward compatibility behaviour
729  if (!theMuonHcalEffects) {
730  savedMuonSimTracks.push_back(myTrack);
731 
732  if (myTrack.onHcal() || myTrack.onVFcal())
733  reconstructHCAL(myTrack, random);
734 
735  return;
736  }
737 
738  if (debug_)
739  LogInfo("FastCalorimetry") << "CalorimetryManager::MuonMipSimulation - track param." << std::endl
740  << " eta = " << moment.eta() << std::endl
741  << " phi = " << moment.phi() << std::endl
742  << " et = " << moment.Et() << std::endl;
743 
744  XYZTLorentzVector trackPosition;
745  if (myTrack.onEcal()) {
746  trackPosition = myTrack.ecalEntrance().vertex();
747  myPart = myTrack.ecalEntrance();
748  } else if (myTrack.onVFcal()) {
749  trackPosition = myTrack.vfcalEntrance().vertex();
750  myPart = myTrack.vfcalEntrance();
751  } else {
752  LogInfo("FastCalorimetry") << " The particle is not in the acceptance " << std::endl;
753  return;
754  }
755 
756  // int onHCAL = hit + 1; - specially for myCalorimeter->hcalProperties(onHCAL)
757  // (below) to get VFcal properties ...
758  // not needed ?
759  // int onHCAL = hit + 1;
760  int onECAL = myTrack.onEcal();
761 
762  //===========================================================================
763 
764  // ECAL and HCAL properties to get
765 
766  //Making ECAL Grid (and segments calculation)
767  XYZPoint caloentrance;
768  XYZVector direction;
769  if (myTrack.onEcal()) {
770  caloentrance = myTrack.ecalEntrance().vertex().Vect();
771  direction = myTrack.ecalEntrance().Vect().Unit();
772  } else if (myTrack.onHcal()) {
773  caloentrance = myTrack.hcalEntrance().vertex().Vect();
774  direction = myTrack.hcalEntrance().Vect().Unit();
775  } else {
776  caloentrance = myTrack.vfcalEntrance().vertex().Vect();
777  direction = myTrack.vfcalEntrance().Vect().Unit();
778  }
779 
780  DetId pivot;
781  if (myTrack.onEcal()) {
782  pivot = myCalorimeter_->getClosestCell(caloentrance, true, myTrack.onEcal() == 1);
783  } else if (myTrack.onHcal()) {
784  pivot = myCalorimeter_->getClosestCell(caloentrance, false, false);
785  }
786 
787  EcalHitMaker myGrid(myCalorimeter_, caloentrance, pivot, pivot.null() ? 0 : myTrack.onEcal(), hdGridSize_, 0, random);
788  // 0 =EM shower -> Unit = X0
789 
790  myGrid.setTrackParameters(direction, 0, myTrack);
791 
792  // Now get the path in the Preshower, ECAL and HCAL along a straight line extrapolation
793  // but only those in the ECAL are used
794 
795  const std::vector<CaloSegment>& segments = myGrid.getSegments();
796  unsigned nsegments = segments.size();
797 
798  int ifirstHcal = -1;
799  int ilastEcal = -1;
800 
802  // // Muon brem in ECAL
803  // MuonBremsstrahlungSimulator* muonBremECAL = 0;
804  // if (theMuonEcalEffects) muonBremECAL = theMuonEcalEffects->muonBremsstrahlungSimulator();
805 
806  for (unsigned iseg = 0; iseg < nsegments && ifirstHcal < 0; ++iseg) {
807  // in the ECAL, there are two types of segments: PbWO4 and GAP
808  float segmentSizeinX0 = segments[iseg].X0length();
809 
810  // Martijn - insert your computations here
811  float energy = 0.0;
812  if (segmentSizeinX0 > 0.001 && segments[iseg].material() == CaloSegment::PbWO4) {
813  // The energy loss simulator
814  float charge = (float)(myTrack.charge());
815  RawParticle p = rawparticle::makeMuon(charge < 0, moment, trackPosition);
816  ParticlePropagator theMuon(p, nullptr, nullptr, mySimEvent->theTable());
817  if (energyLossECAL) {
818  energyLossECAL->updateState(theMuon, segmentSizeinX0, random);
819  energy = energyLossECAL->deltaMom().E();
820  moment -= energyLossECAL->deltaMom();
821  }
822  }
823  // that's all for ECAL, Florian
824  // Save the hit only if it is a crystal
825  if (segments[iseg].material() == CaloSegment::PbWO4) {
826  myGrid.getPads(segments[iseg].sX0Entrance() + segmentSizeinX0 * 0.5);
827  myGrid.setSpotEnergy(energy);
828  myGrid.addHit(0., 0.);
829  ilastEcal = iseg;
830  }
831  // Check for end of loop:
832  if (segments[iseg].material() == CaloSegment::HCAL) {
833  ifirstHcal = iseg;
834  }
835  }
836 
837  // Build the FAMOS HCAL
838  HcalHitMaker myHcalHitMaker(myGrid, (unsigned)2);
839  // float mipenergy=0.1;
840  // Create the helix with the stepping helix propagator
841  // to add a hit, just do
842  // myHcalHitMaker.setSpotEnergy(mipenergy);
843  // math::XYZVector hcalEntrance;
844  // if(ifirstHcal>=0) hcalEntrance=segments[ifirstHcal].entrance();
845  // myHcalHitMaker.addHit(hcalEntrance);
849  int ilastHcal = -1;
850  float mipenergy = 0.0;
851 
853  // // Muon Brem effect
854  // MuonBremsstrahlungSimulator* muonBremHCAL = 0;
855  // if (theMuonHcalEffects) muonBremHCAL = theMuonHcalEffects->muonBremsstrahlungSimulator();
856 
857  if (ifirstHcal > 0 && energyLossHCAL) {
858  for (unsigned iseg = ifirstHcal; iseg < nsegments; ++iseg) {
859  float segmentSizeinX0 = segments[iseg].X0length();
860  if (segments[iseg].material() == CaloSegment::HCAL) {
861  ilastHcal = iseg;
862  if (segmentSizeinX0 > 0.001) {
863  // The energy loss simulator
864  float charge = (float)(myTrack.charge());
865  RawParticle p = rawparticle::makeMuon(charge < 0, moment, trackPosition);
866  ParticlePropagator theMuon(p, nullptr, nullptr, mySimEvent->theTable());
867  energyLossHCAL->updateState(theMuon, segmentSizeinX0, random);
868  mipenergy = energyLossHCAL->deltaMom().E();
869  moment -= energyLossHCAL->deltaMom();
870  myHcalHitMaker.setSpotEnergy(mipenergy);
871  myHcalHitMaker.addHit(segments[iseg].entrance());
872  }
873  }
874  }
875  }
876 
877  // Copy the muon SimTrack (Only for Energy loss)
878  FSimTrack muonTrack(myTrack);
879  if (energyLossHCAL && ilastHcal >= 0) {
880  math::XYZVector hcalExit = segments[ilastHcal].exit();
881  muonTrack.setTkPosition(hcalExit);
882  muonTrack.setTkMomentum(moment);
883  } else if (energyLossECAL && ilastEcal >= 0) {
884  math::XYZVector ecalExit = segments[ilastEcal].exit();
885  muonTrack.setTkPosition(ecalExit);
886  muonTrack.setTkMomentum(moment);
887  } // else just leave tracker surface position and momentum...
888 
889  muonSimTracks.push_back(muonTrack);
890 
891  // no need to change below this line
892  std::map<CaloHitID, float>::const_iterator mapitr;
893  std::map<CaloHitID, float>::const_iterator endmapitr;
894  if (myTrack.onEcal() > 0) {
895  // Save ECAL hits
896  updateECAL(myGrid.getHits(), onECAL, myTrack.id());
897  }
898 
899  // Save HCAL hits
900  updateHCAL(myHcalHitMaker.getHits(), myTrack.id());
901 
902  if (debug_)
903  LogInfo("FastCalorimetry") << std::endl << " FASTEnergyReconstructor::MipShowerSimulation finished " << std::endl;
904 }
905 
907  edm::ParameterSet ECALparameters = fastCalo.getParameter<edm::ParameterSet>("ECAL");
908 
909  evtsToDebug_ = fastCalo.getUntrackedParameter<std::vector<unsigned int> >("EvtsToDebug", std::vector<unsigned>());
910  debug_ = fastCalo.getUntrackedParameter<bool>("Debug");
911 
912  bFixedLength_ = ECALparameters.getParameter<bool>("bFixedLength");
913 
914  gridSize_ = ECALparameters.getParameter<int>("GridSize");
915  spotFraction_ = ECALparameters.getParameter<double>("SpotFraction");
916  pulledPadSurvivalProbability_ = ECALparameters.getParameter<double>("FrontLeakageProbability");
917  crackPadSurvivalProbability_ = ECALparameters.getParameter<double>("GapLossProbability");
918  theCoreIntervals_ = ECALparameters.getParameter<std::vector<double> >("CoreIntervals");
919  theTailIntervals_ = ECALparameters.getParameter<std::vector<double> >("TailIntervals");
920 
921  RCFactor_ = ECALparameters.getParameter<double>("RCFactor");
922  RTFactor_ = ECALparameters.getParameter<double>("RTFactor");
923  //changed after tuning - Feb-July - Shilpi Jain
924  // radiusFactor_ = ECALparameters.getParameter<double>("RadiusFactor");
925  radiusFactorEE_ = ECALparameters.getParameter<double>("RadiusFactorEE");
926  radiusFactorEB_ = ECALparameters.getParameter<double>("RadiusFactorEB");
927  //(end of) changed after tuning - Feb-July - Shilpi Jain
928  radiusPreshowerCorrections_ = ECALparameters.getParameter<std::vector<double> >("RadiusPreshowerCorrections");
929  aTerm = 1. + radiusPreshowerCorrections_[1] * radiusPreshowerCorrections_[0];
930  bTerm = radiusPreshowerCorrections_[0];
931  mipValues_ = ECALparameters.getParameter<std::vector<double> >("MipsinGeV");
932  simulatePreshower_ = ECALparameters.getParameter<bool>("SimulatePreshower");
933 
934  if (gridSize_ < 1)
935  gridSize_ = 7;
936  if (pulledPadSurvivalProbability_ < 0. || pulledPadSurvivalProbability_ > 1)
937  pulledPadSurvivalProbability_ = 1.;
938  if (crackPadSurvivalProbability_ < 0. || crackPadSurvivalProbability_ > 1)
939  crackPadSurvivalProbability_ = 0.9;
940 
941  LogInfo("FastCalorimetry") << " Fast ECAL simulation parameters " << std::endl;
942  LogInfo("FastCalorimetry") << " =============================== " << std::endl;
943  if (simulatePreshower_)
944  LogInfo("FastCalorimetry") << " The preshower is present " << std::endl;
945  else
946  LogInfo("FastCalorimetry") << " The preshower is NOT present " << std::endl;
947  LogInfo("FastCalorimetry") << " Grid Size : " << gridSize_ << std::endl;
948  if (spotFraction_ > 0.)
949  LogInfo("FastCalorimetry") << " Spot Fraction : " << spotFraction_ << std::endl;
950  else {
951  LogInfo("FastCalorimetry") << " Core of the shower " << std::endl;
952  for (unsigned ir = 0; ir < theCoreIntervals_.size() / 2; ++ir) {
953  LogInfo("FastCalorimetry") << " r < " << theCoreIntervals_[ir * 2] << " R_M : " << theCoreIntervals_[ir * 2 + 1]
954  << " ";
955  }
956  LogInfo("FastCalorimetry") << std::endl;
957 
958  LogInfo("FastCalorimetry") << " Tail of the shower " << std::endl;
959  for (unsigned ir = 0; ir < theTailIntervals_.size() / 2; ++ir) {
960  LogInfo("FastCalorimetry") << " r < " << theTailIntervals_[ir * 2] << " R_M : " << theTailIntervals_[ir * 2 + 1]
961  << " ";
962  }
963  //changed after tuning - Feb-July - Shilpi Jain
964  LogInfo("FastCalorimetry") << "Radius correction factors: EB & EE " << radiusFactorEB_ << " : " << radiusFactorEE_
965  << std::endl;
966  //(end of) changed after tuning - Feb-July - Shilpi Jain
967  LogInfo("FastCalorimetry") << std::endl;
968  if (mipValues_.size() > 2) {
969  LogInfo("FastCalorimetry") << "Improper number of parameters for the preshower ; using 95keV" << std::endl;
970  mipValues_.clear();
971  mipValues_.resize(2, 0.000095);
972  }
973  }
974 
975  LogInfo("FastCalorimetry") << " FrontLeakageProbability : " << pulledPadSurvivalProbability_ << std::endl;
976  LogInfo("FastCalorimetry") << " GapLossProbability : " << crackPadSurvivalProbability_ << std::endl;
977 
978  // RespCorrP: p (momentum), ECAL and HCAL corrections = f(p)
979  edm::ParameterSet CalorimeterParam = fastCalo.getParameter<edm::ParameterSet>("CalorimeterProperties");
980 
981  rsp = CalorimeterParam.getParameter<std::vector<double> >("RespCorrP");
982  LogInfo("FastCalorimetry") << " RespCorrP (rsp) size " << rsp.size() << std::endl;
983 
984  if (rsp.size() % 3 != 0) {
985  LogInfo("FastCalorimetry") << " RespCorrP size is wrong -> no corrections applied !!!" << std::endl;
986 
987  p_knots.push_back(14000.);
988  k_e.push_back(1.);
989  k_h.push_back(1.);
990  } else {
991  for (unsigned i = 0; i < rsp.size(); i += 3) {
992  LogInfo("FastCalorimetry") << "i = " << i / 3 << " p = " << rsp[i] << " k_e(p) = " << rsp[i + 1]
993  << " k_e(p) = " << rsp[i + 2] << std::endl;
994 
995  p_knots.push_back(rsp[i]);
996  k_e.push_back(rsp[i + 1]);
997  k_h.push_back(rsp[i + 2]);
998  }
999  }
1000 
1001  //FR
1002  edm::ParameterSet HCALparameters = fastCalo.getParameter<edm::ParameterSet>("HCAL");
1003  optionHDSim_ = HCALparameters.getParameter<int>("SimOption");
1004  hdGridSize_ = HCALparameters.getParameter<int>("GridSize");
1005  hdSimMethod_ = HCALparameters.getParameter<int>("SimMethod");
1006  //RF
1007 
1008  EcalDigitizer_ = ECALparameters.getUntrackedParameter<bool>("Digitizer", false);
1009  HcalDigitizer_ = HCALparameters.getUntrackedParameter<bool>("Digitizer", false);
1010  samplingHBHE_ = HCALparameters.getParameter<std::vector<double> >("samplingHBHE");
1011  samplingHF_ = HCALparameters.getParameter<std::vector<double> >("samplingHF");
1012  samplingHO_ = HCALparameters.getParameter<std::vector<double> >("samplingHO");
1013  ietaShiftHB_ = HCALparameters.getParameter<int>("ietaShiftHB");
1014  ietaShiftHE_ = HCALparameters.getParameter<int>("ietaShiftHE");
1015  ietaShiftHF_ = HCALparameters.getParameter<int>("ietaShiftHF");
1016  ietaShiftHO_ = HCALparameters.getParameter<int>("ietaShiftHO");
1017  timeShiftHB_ = HCALparameters.getParameter<std::vector<double> >("timeShiftHB");
1018  timeShiftHE_ = HCALparameters.getParameter<std::vector<double> >("timeShiftHE");
1019  timeShiftHF_ = HCALparameters.getParameter<std::vector<double> >("timeShiftHF");
1020  timeShiftHO_ = HCALparameters.getParameter<std::vector<double> >("timeShiftHO");
1021 
1022  // FastHFShowerLibrary
1023  edm::ParameterSet m_HS = fastCalo.getParameter<edm::ParameterSet>("HFShowerLibrary");
1024  useShowerLibrary = m_HS.getUntrackedParameter<bool>("useShowerLibrary", false);
1025  useCorrectionSL = m_HS.getUntrackedParameter<bool>("useCorrectionSL", false);
1026 }
1027 
1029  int sizeP = p_knots.size();
1030 
1031  if (sizeP <= 1) {
1032  ecorr = 1.;
1033  hcorr = 1.;
1034  } else {
1035  int ip = -1;
1036  for (int i = 0; i < sizeP; i++) {
1037  if (p < p_knots[i]) {
1038  ip = i;
1039  break;
1040  }
1041  }
1042  if (ip == 0) {
1043  ecorr = k_e[0];
1044  hcorr = k_h[0];
1045  } else {
1046  if (ip == -1) {
1047  ecorr = k_e[sizeP - 1];
1048  hcorr = k_h[sizeP - 1];
1049  } else {
1050  double x1 = p_knots[ip - 1];
1051  double x2 = p_knots[ip];
1052  double y1 = k_e[ip - 1];
1053  double y2 = k_e[ip];
1054 
1055  ecorr = (y1 + (y2 - y1) * (p - x1) / (x2 - x1));
1056 
1057  y1 = k_h[ip - 1];
1058  y2 = k_h[ip];
1059  hcorr = (y1 + (y2 - y1) * (p - x1) / (x2 - x1));
1060  }
1061  }
1062  }
1063 
1064  if (debug_)
1065  LogInfo("FastCalorimetry") << " p, ecorr, hcorr = " << p << " " << ecorr << " " << hcorr << std::endl;
1066 }
1067 
1068 void CalorimetryManager::updateECAL(const std::map<CaloHitID, float>& hitMap, int onEcal, int trackID, float corr) {
1069  std::map<CaloHitID, float>::const_iterator mapitr;
1070  std::map<CaloHitID, float>::const_iterator endmapitr = hitMap.end();
1071  if (onEcal == 1) {
1072  EBMapping_.reserve(EBMapping_.size() + hitMap.size());
1073  endmapitr = hitMap.end();
1074  for (mapitr = hitMap.begin(); mapitr != endmapitr; ++mapitr) {
1075  //correct energy
1076  float energy = mapitr->second;
1077  energy *= corr;
1078 
1079  //make finalized CaloHitID
1080  CaloHitID current_id(mapitr->first.unitID(), mapitr->first.timeSlice(), trackID);
1081 
1082  EBMapping_.push_back(std::pair<CaloHitID, float>(current_id, energy));
1083  }
1084  } else if (onEcal == 2) {
1085  EEMapping_.reserve(EEMapping_.size() + hitMap.size());
1086  endmapitr = hitMap.end();
1087  for (mapitr = hitMap.begin(); mapitr != endmapitr; ++mapitr) {
1088  //correct energy
1089  float energy = mapitr->second;
1090  energy *= corr;
1091 
1092  //make finalized CaloHitID
1093  CaloHitID current_id(mapitr->first.unitID(), mapitr->first.timeSlice(), trackID);
1094 
1095  EEMapping_.push_back(std::pair<CaloHitID, float>(current_id, energy));
1096  }
1097  }
1098 }
1099 
1100 void CalorimetryManager::updateHCAL(const std::map<CaloHitID, float>& hitMap, int trackID, float corr) {
1101  std::vector<double> hfcorrEm = myHDResponse_->getCorrHFem();
1102  std::vector<double> hfcorrHad = myHDResponse_->getCorrHFhad();
1103  std::map<CaloHitID, float>::const_iterator mapitr;
1104  std::map<CaloHitID, float>::const_iterator endmapitr = hitMap.end();
1105  HMapping_.reserve(HMapping_.size() + hitMap.size());
1106  for (mapitr = hitMap.begin(); mapitr != endmapitr; ++mapitr) {
1107  //correct energy
1108  float energy = mapitr->second;
1109  energy *= corr;
1110 
1111  float time = mapitr->first.timeSlice();
1112  //put energy into uncalibrated state for digitizer && correct timing
1113  if (HcalDigitizer_) {
1114  HcalDetId hdetid = HcalDetId(mapitr->first.unitID());
1115  if (hdetid.subdetId() == HcalBarrel) {
1116  energy /= samplingHBHE_[hdetid.ietaAbs() - 1]; //re-convert to GeV
1117  time = timeShiftHB_[hdetid.ietaAbs() - ietaShiftHB_];
1118  } else if (hdetid.subdetId() == HcalEndcap) {
1119  energy /= samplingHBHE_[hdetid.ietaAbs() - 1]; //re-convert to GeV
1120  time = timeShiftHE_[hdetid.ietaAbs() - ietaShiftHE_];
1121  } else if (hdetid.subdetId() == HcalForward) {
1122  if (useShowerLibrary) {
1123  if (useCorrectionSL) {
1124  if (hdetid.depth() == 1 or hdetid.depth() == 3)
1125  energy *= hfcorrEm[hdetid.ietaAbs() - ietaShiftHF_];
1126  if (hdetid.depth() == 2 or hdetid.depth() == 4)
1127  energy *= hfcorrHad[hdetid.ietaAbs() - ietaShiftHF_];
1128  }
1129  } else {
1130  if (hdetid.depth() == 1 or hdetid.depth() == 3)
1131  energy *= samplingHF_[0];
1132  if (hdetid.depth() == 2 or hdetid.depth() == 4)
1133  energy *= samplingHF_[1];
1134  time = timeShiftHF_[hdetid.ietaAbs() - ietaShiftHF_];
1135  }
1136  } else if (hdetid.subdetId() == HcalOuter) {
1137  energy /= samplingHO_[hdetid.ietaAbs() - 1];
1138  time = timeShiftHO_[hdetid.ietaAbs() - ietaShiftHO_];
1139  }
1140  }
1141 
1142  //make finalized CaloHitID
1143  CaloHitID current_id(mapitr->first.unitID(), time, trackID);
1144  HMapping_.push_back(std::pair<CaloHitID, float>(current_id, energy));
1145  }
1146 }
1147 
1148 void CalorimetryManager::updatePreshower(const std::map<CaloHitID, float>& hitMap, int trackID, float corr) {
1149  std::map<CaloHitID, float>::const_iterator mapitr;
1150  std::map<CaloHitID, float>::const_iterator endmapitr = hitMap.end();
1151  ESMapping_.reserve(ESMapping_.size() + hitMap.size());
1152  for (mapitr = hitMap.begin(); mapitr != endmapitr; ++mapitr) {
1153  //correct energy
1154  float energy = mapitr->second;
1155  energy *= corr;
1156 
1157  //make finalized CaloHitID
1158  CaloHitID current_id(mapitr->first.unitID(), mapitr->first.timeSlice(), trackID);
1159 
1160  ESMapping_.push_back(std::pair<CaloHitID, float>(current_id, energy));
1161  }
1162 }
1163 
1165  c.reserve(c.size() + EBMapping_.size());
1166  for (unsigned i = 0; i < EBMapping_.size(); i++) {
1167  c.push_back(PCaloHit(EBDetId::unhashIndex(EBMapping_[i].first.unitID()),
1168  EBMapping_[i].second,
1169  EBMapping_[i].first.timeSlice(),
1170  EBMapping_[i].first.trackID()));
1171  }
1172 }
1173 
1175  c.reserve(c.size() + EEMapping_.size());
1176  for (unsigned i = 0; i < EEMapping_.size(); i++) {
1177  c.push_back(PCaloHit(EEDetId::unhashIndex(EEMapping_[i].first.unitID()),
1178  EEMapping_[i].second,
1179  EEMapping_[i].first.timeSlice(),
1180  EEMapping_[i].first.trackID()));
1181  }
1182 }
1183 
1185  c.reserve(c.size() + HMapping_.size());
1186  for (unsigned i = 0; i < HMapping_.size(); i++) {
1187  c.push_back(PCaloHit(DetId(HMapping_[i].first.unitID()),
1188  HMapping_[i].second,
1189  HMapping_[i].first.timeSlice(),
1190  HMapping_[i].first.trackID()));
1191  }
1192 }
1193 
1195  c.reserve(c.size() + ESMapping_.size());
1196  for (unsigned i = 0; i < ESMapping_.size(); i++) {
1197  c.push_back(PCaloHit(ESMapping_[i].first.unitID(),
1198  ESMapping_[i].second,
1199  ESMapping_[i].first.timeSlice(),
1200  ESMapping_[i].first.trackID()));
1201  }
1202 }
1203 
1204 // The main danger in this method is to screw up to relationships between particles
1205 // So, the muon FSimTracks created by FSimEvent.cc are simply to be updated
1207  unsigned size = muons.size();
1208  for (unsigned i = 0; i < size; ++i) {
1209  int id = muons[i].trackId();
1210  if (!(abs(muons[i].type()) == 13 || abs(muons[i].type()) == 1000024 ||
1211  (abs(muons[i].type()) > 1000100 && abs(muons[i].type()) < 1999999)))
1212  continue;
1213  // identify the corresponding muon in the local collection
1214 
1215  std::vector<FSimTrack>::const_iterator itcheck =
1216  find_if(muonSimTracks.begin(), muonSimTracks.end(), FSimTrackEqual(id));
1217  if (itcheck != muonSimTracks.end()) {
1218  muons[i].setTkPosition(itcheck->trackerSurfacePosition());
1219  muons[i].setTkMomentum(itcheck->trackerSurfaceMomentum());
1220  }
1221  }
1222 }
1223 
1225  c.reserve(int(0.2 * muonSimTracks.size() + 0.2 * savedMuonSimTracks.size() + 0.5));
1226  for (const auto& track : muonSimTracks) {
1227  if (track.momentum().perp2() > 1.0 && fabs(track.momentum().eta()) < 3.0 && track.isGlobal())
1228  c.push_back(track);
1229  }
1230  for (const auto& track : savedMuonSimTracks) {
1231  if (track.momentum().perp2() > 1.0 && fabs(track.momentum().eta()) < 3.0 && track.isGlobal())
1232  c.push_back(track);
1233  }
1234  c.shrink_to_fit();
1235 }
size
Write out results.
type
Definition: HCALResponse.h:21
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:40
T getUntrackedParameter(std::string const &, T const &) const
std::vector< double > rsp
std::vector< double > k_h
bool noEndVertex() const
no end vertex
RawParticle myElec
A few pointers to save time.
float charge() const
charge
Definition: FSimTrack.h:56
void reconstructHCAL(const FSimTrack &myTrack, RandomEngineAndDistribution const *)
std::vector< PCaloHit > PCaloHitContainer
XYZVector Vect() const
the momentum threevector
Definition: RawParticle.h:323
const double GeV
Definition: MathUtil.h:16
const ECALProperties * ecalProperties(int onEcal) const
ECAL properties.
Definition: Calorimeter.cc:78
double responseHCAL(int _mip, double energy, double eta, int partype, RandomEngineAndDistribution const *)
const RawParticle & vfcalEntrance() const
The particle at VFCAL entrance.
Definition: FSimTrack.h:149
void updateHCAL(const std::map< CaloHitID, float > &hitMap, int trackID=0, float corr=1.0)
int onLayer2() const
Definition: FSimTrack.h:106
GflashPiKShowerProfile * thePiKProfile
double flatShoot(double xmin=0.0, double xmax=1.0) const
std::vector< double > timeShiftHO_
void updatePreshower(const std::map< CaloHitID, float > &hitMap, int trackID=0, float corr=1.0)
GflashTrajectory * getHelix()
Definition: GflashShowino.h:33
constexpr bool null() const
is this a null id ?
Definition: DetId.h:59
MaterialEffects * theMuonEcalEffects
MaterialEffects * theMuonHcalEffects
void recoHFShowerLibrary(const FSimTrack &myTrack)
const XYZTLorentzVector & momentum() const
Temporary (until move of SimTrack to Mathcore) - No! Actually very useful.
Definition: FSimTrack.h:209
void setCrackPadSurvivalProbability(double val)
Definition: EcalHitMaker.h:129
double getHCALEnergyResponse(double e, int hit, RandomEngineAndDistribution const *)
GflashHadronShowerProfile * theProfile
const HepPDT::ParticleDataTable * theTable() const
Get the pointer to the particle data table.
Definition: FBaseSimEvent.h:54
#define nullptr
const RawParticle & layer1Entrance() const
The particle at Preshower Layer 1.
Definition: FSimTrack.h:137
const std::map< CaloHitID, float > & getHits() override
bool compute()
Compute the shower longitudinal and lateral development.
Definition: HDShower.cc:453
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
bool exists(std::string const &parameterName) const
checks if a parameter exists
void correctHF(double e, int type)
void reconstructTrack(FSimTrack &myTrack, RandomEngineAndDistribution const *)
void setPreshower(PreshowerHitMaker *const myPresh)
set the preshower address
Definition: EMShower.cc:654
void loadFromPreshower(edm::PCaloHitContainer &c) const
Definition: weight.py:1
void updateState(ParticlePropagator &myTrack, double radlen, RandomEngineAndDistribution const *)
Compute the material effect (calls the sub class)
void updateECAL(const std::map< CaloHitID, float > &hitMap, int onEcal, int trackID=0, float corr=1.0)
const CaloSubdetectorGeometry * getHcalGeometry() const
Definition: Calorimeter.h:56
void MuonMipSimulation(const FSimTrack &myTrack, RandomEngineAndDistribution const *)
void readParameters(const edm::ParameterSet &fastCalo)
RawParticle makeMuon(bool isParticle, const math::XYZTLorentzVector &p, const math::XYZTLorentzVector &xStart)
Definition: makeMuon.cc:20
const PreshowerLayer1Properties * layer1Properties(int onLayer1) const
Preshower Layer1 properties.
Definition: Calorimeter.cc:103
static EEDetId unhashIndex(int hi)
Definition: EEDetId.cc:65
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
static double constexpr eMass
Electron mass[GeV].
Definition: Constants.h:13
std::vector< double > timeShiftHF_
std::vector< double > p_knots
void compute()
Compute the shower longitudinal and lateral development.
Definition: EMShower.cc:263
int onEcal() const
Definition: FSimTrack.h:111
std::vector< FSimTrack > muonSimTracks
double radLenIncm() const override
Radiation length in cm.
U second(std::pair< T, U > const &p)
void HDShowerSimulation(const FSimTrack &myTrack, RandomEngineAndDistribution const *)
Hadronic Shower Simulation.
const LandauFluctuationGenerator * aLandauGenerator
std::vector< double > samplingHF_
bool compute()
Compute the shower longitudinal and lateral development.
Definition: HFShower.cc:461
void setRadiusFactor(double r)
Definition: EcalHitMaker.h:125
bool addHit(double r, double phi, unsigned layer=0) override
add the hit in the HCAL in local coordinates
Definition: HcalHitMaker.cc:30
int depth() const
get the tower depth
Definition: HcalDetId.h:164
const XYZTLorentzVector & momentum() const
the momentum fourvector
Definition: RawParticle.h:321
void reconstruct(RandomEngineAndDistribution const *)
CalorimeterNumber getCalorimeterNumber(const Gflash3Vector &position)
double e() const
energy of the momentum
Definition: RawParticle.h:305
GflashAntiProtonShowerProfile * theAntiProtonProfile
math::XYZVector XYZVector
int getmip()
Definition: HDShower.h:39
GflashProtonShowerProfile * theProtonProfile
const PreshowerLayer2Properties * layer2Properties(int onLayer2) const
Preshower Layer2 properties.
Definition: Calorimeter.cc:110
double getPathLengthAtShower()
Definition: GflashShowino.h:29
const std::map< CaloHitID, float > & getHits() override
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
void setTrackParameters(const XYZNormal &normal, double X0depthoffset, const FSimTrack &theTrack)
std::vector< std::pair< CaloHitID, float > > ESMapping_
void loadFromEcalEndcap(edm::PCaloHitContainer &c) const
const XYZTLorentzVector & deltaMom() const
Returns the actual energy lost.
std::vector< double > radiusPreshowerCorrections_
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
const HCALProperties * hcalProperties(int onHcal) const
HCAL properties.
Definition: Calorimeter.cc:88
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double getMIPfraction(double energy, double eta)
vec1 & getCorrHFhad()
Definition: HCALResponse.h:48
const RawParticle & ecalEntrance() const
The particle at ECAL entrance.
Definition: FSimTrack.h:143
int onVFcal() const
Definition: FSimTrack.h:121
void setSpotEnergy(double e) override
Set the spot energy.
Definition: HcalHitMaker.h:26
const std::map< CaloHitID, float > & getHitsMap()
std::vector< double > theTailIntervals_
std::vector< double > samplingHBHE_
void loadMuonSimTracks(edm::SimTrackContainer &m) const
unsigned int nTracks() const
Number of tracks.
Definition: FSimEvent.cc:24
JetCorrectorParameters corr
Definition: classes.h:5
HCALResponse * myHDResponse_
vec1 & getCorrHFem()
Definition: HCALResponse.h:47
CaloGeometryHelper * myCalorimeter_
std::vector< double > timeShiftHB_
int ietaAbs() const
get the absolute value of the cell ieta
Definition: HcalDetId.h:148
std::vector< std::pair< CaloHitID, float > > HMapping_
void initialize(RandomEngineAndDistribution const *random)
const double intLength[kNumberCalorimeter]
std::vector< double > theCoreIntervals_
void getGflashTrajectoryPoint(GflashTrajectoryPoint &point, double s) const
edm::EventID id() const
Method to return the EventId.
Definition: FSimEvent.cc:20
const XYZTLorentzVector & vertex() const
the vertex fourvector
Definition: RawParticle.h:320
Definition: DetId.h:17
void loadFromHcal(edm::PCaloHitContainer &c) const
GammaFunctionGenerator * aGammaGenerator
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:31
CLHEP::Hep3Vector Gflash3Vector
Definition: Gflash3Vector.h:6
std::vector< double > mipValues_
static std::vector< std::pair< int, float > > myZero_
std::vector< FSimTrack > savedMuonSimTracks
void setTkPosition(const math::XYZVectorD &pos)
Definition: SimTrack.h:41
std::vector< std::pair< CaloHitID, float > > EEMapping_
double getPathLengthOnEcal()
Definition: GflashShowino.h:28
void SetRandom(const RandomEngineAndDistribution *)
static EBDetId unhashIndex(int hi)
get a DetId from a compact index for arrays
Definition: EBDetId.h:110
math::XYZVector XYZPoint
int type() const
particle type (HEP PDT convension)
Definition: CoreSimTrack.h:22
HLT enums.
void harvestMuonSimTracks(edm::SimTrackContainer &m) const
std::vector< unsigned int > evtsToDebug_
int id() const
the index in FBaseSimEvent and other vectors
Definition: FSimTrack.h:96
bool preshowerPresent() const
void setPulledPadSurvivalProbability(double val)
Definition: EcalHitMaker.h:127
const RawParticle & layer2Entrance() const
The particle at Preshower Layer 2.
Definition: FSimTrack.h:140
int onLayer1() const
Definition: FSimTrack.h:101
std::vector< double > k_e
std::vector< GflashHit > & getGflashHitList()
int onHcal() const
Definition: FSimTrack.h:116
void setTkMomentum(const math::XYZTLorentzVectorD &mom)
Definition: SimTrack.h:43
double eta() const
Definition: RawParticle.h:279
const RawParticle & hcalEntrance() const
The particle at HCAL entrance.
Definition: FSimTrack.h:146
std::vector< std::pair< CaloHitID, float > > EBMapping_
std::vector< double > samplingHO_
void setPreshowerPresent(bool ps)
Definition: EcalHitMaker.h:132
FastHFShowerLibrary * theHFShowerLibrary
DetId getClosestCell(const XYZPoint &point, bool ecal, bool central) const
math::XYZVector XYZVector
Definition: RawParticle.h:26
EnergyLossSimulator * energyLossSimulator() const
Return the Energy Loss engine.
std::unique_ptr< KKCorrectionFactors > ecalCorrection
HSParameters * myHSParameters_
std::vector< SimTrack > SimTrackContainer
void setVertex(const XYZTLorentzVector &vtx)
set the vertex
Definition: RawParticle.h:325
bool setDepth(double, bool inCm=false)
set the depth in X0 or Lambda0 units depending on showerType
void initialize(int showerType, double energy, double globalTime, double charge, Gflash3Vector &position, Gflash3Vector &momentum)
bool computeShower()
Definition: HDRShower.cc:41
math::XYZTLorentzVector XYZTLorentzVector
Definition: RawParticle.h:25
void print() const
print the FBaseSimEvent in an intelligible way
void setHcal(HcalHitMaker *const myHcal)
set the HCAL address
Definition: EMShower.cc:661
void setGrid(EcalHitMaker *const myGrid)
set the grid address
Definition: EMShower.h:58
void EMShowerSimulation(const FSimTrack &myTrack, RandomEngineAndDistribution const *)
const std::map< CaloHitID, float > & getHits() override
Definition: HcalHitMaker.h:35
std::vector< double > timeShiftHE_
FSimTrack & track(int id) const
Return track with given Id.
double getMaximumOfShower() const
get the depth of the centre of gravity of the shower(s)
Definition: EMShower.h:55
void loadFromEcalBarrel(edm::PCaloHitContainer &c) const