Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:11:14 

0001 //updated by Reza Goldouzian
0002 //Framework headers
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/Utilities/interface/EDMException.h"
0006 
0007 // Fast Simulation headers
0008 #include "FastSimulation/ParticlePropagator/interface/ParticlePropagator.h"
0009 #include "FastSimulation/Calorimetry/interface/CalorimetryManager.h"
0010 #include "FastSimulation/Event/interface/FSimEvent.h"
0011 #include "FastSimulation/Event/interface/FSimTrack.h"
0012 #include "FastSimulation/ShowerDevelopment/interface/EMECALShowerParametrization.h"
0013 #include "FastSimulation/ShowerDevelopment/interface/EMShower.h"
0014 #include "FastSimulation/ShowerDevelopment/interface/HDShowerParametrization.h"
0015 #include "FastSimulation/ShowerDevelopment/interface/HDShower.h"
0016 #include "FastSimulation/ShowerDevelopment/interface/HFShower.h"
0017 #include "FastSimulation/ShowerDevelopment/interface/HDRShower.h"
0018 #include "FastSimulation/ShowerDevelopment/interface/HSParameters.h"
0019 #include "FastSimulation/CaloGeometryTools/interface/CaloGeometryHelper.h"
0020 //#include "FastSimulation/Utilities/interface/Histos.h"
0021 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0022 #include "FastSimulation/Utilities/interface/GammaFunctionGenerator.h"
0023 #include "FastSimulation/Utilities/interface/LandauFluctuationGenerator.h"
0024 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0025 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0026 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0027 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0028 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0029 #include "FastSimulation/Event/interface/FSimTrackEqual.h"
0030 // New headers for Muon Mip Simulation
0031 #include "FastSimulation/MaterialEffects/interface/MaterialEffects.h"
0032 #include "FastSimulation/MaterialEffects/interface/EnergyLossSimulator.h"
0033 // Muon Brem
0034 #include "FastSimulation/MaterialEffects/interface/MuonBremsstrahlungSimulator.h"
0035 
0036 //Gflash Hadronic Model
0037 #include "SimGeneral/GFlash/interface/GflashHadronShowerProfile.h"
0038 #include "SimGeneral/GFlash/interface/GflashPiKShowerProfile.h"
0039 #include "SimGeneral/GFlash/interface/GflashProtonShowerProfile.h"
0040 #include "SimGeneral/GFlash/interface/GflashAntiProtonShowerProfile.h"
0041 #include "SimGeneral/GFlash/interface/GflashTrajectoryPoint.h"
0042 #include "SimGeneral/GFlash/interface/GflashHit.h"
0043 #include "SimGeneral/GFlash/interface/Gflash3Vector.h"
0044 
0045 //FastHFShowerLibrary
0046 #include "FastSimulation/ShowerDevelopment/interface/FastHFShowerLibrary.h"
0047 
0048 // STL headers
0049 #include <memory>
0050 
0051 #include <iostream>
0052 #include <vector>
0053 
0054 //CMSSW headers
0055 #include "DataFormats/DetId/interface/DetId.h"
0056 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0057 //#include "DataFormats/EcalDetId/interface/EcalDetId.h"
0058 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0059 
0060 //ROOT headers
0061 #include "TROOT.h"
0062 #include "TH1.h"
0063 
0064 using namespace edm;
0065 
0066 typedef math::XYZVector XYZVector;
0067 typedef math::XYZVector XYZPoint;
0068 
0069 std::vector<std::pair<int, float> > CalorimetryManager::myZero_ =
0070     std::vector<std::pair<int, float> >(1, std::pair<int, float>(0, 0.));
0071 
0072 CalorimetryManager::CalorimetryManager() : myCalorimeter_(nullptr), initialized_(false) { ; }
0073 
0074 CalorimetryManager::CalorimetryManager(FSimEvent* aSimEvent,
0075                                        const edm::ParameterSet& fastCalo,
0076                                        const edm::ParameterSet& fastMuECAL,
0077                                        const edm::ParameterSet& fastMuHCAL,
0078                                        const edm::ParameterSet& parGflash,
0079                                        edm::ConsumesCollector&& iC)
0080     : mySimEvent(aSimEvent),
0081       initialized_(false),
0082       theMuonEcalEffects(nullptr),
0083       theMuonHcalEffects(nullptr),
0084       bFixedLength_(false) {
0085   aLandauGenerator = new LandauFluctuationGenerator;
0086   aGammaGenerator = new GammaFunctionGenerator;
0087 
0088   //Gflash
0089   theProfile = new GflashHadronShowerProfile(parGflash);
0090   thePiKProfile = new GflashPiKShowerProfile(parGflash);
0091   theProtonProfile = new GflashProtonShowerProfile(parGflash);
0092   theAntiProtonProfile = new GflashAntiProtonShowerProfile(parGflash);
0093 
0094   // FastHFShowerLibrary
0095   theHFShowerLibrary = new FastHFShowerLibrary(fastCalo, std::move(iC));
0096 
0097   readParameters(fastCalo);
0098 
0099   myCalorimeter_ = new CaloGeometryHelper(fastCalo);
0100   myHDResponse_ = new HCALResponse(fastCalo.getParameter<edm::ParameterSet>("HCALResponse"));
0101   myHSParameters_ = new HSParameters(fastCalo.getParameter<edm::ParameterSet>("HSParameters"));
0102 
0103   // Material Effects for Muons in ECAL (only EnergyLoss implemented so far)
0104   if (fastMuECAL.getParameter<bool>("PairProduction") || fastMuECAL.getParameter<bool>("Bremsstrahlung") ||
0105       fastMuECAL.getParameter<bool>("MuonBremsstrahlung") || fastMuECAL.getParameter<bool>("EnergyLoss") ||
0106       fastMuECAL.getParameter<bool>("MultipleScattering"))
0107     theMuonEcalEffects = new MaterialEffects(fastMuECAL);
0108 
0109   // Material Effects for Muons in HCAL (only EnergyLoss implemented so far)
0110   if (fastMuHCAL.getParameter<bool>("PairProduction") || fastMuHCAL.getParameter<bool>("Bremsstrahlung") ||
0111       fastMuHCAL.getParameter<bool>("MuonBremsstrahlung") || fastMuHCAL.getParameter<bool>("EnergyLoss") ||
0112       fastMuHCAL.getParameter<bool>("MultipleScattering"))
0113     theMuonHcalEffects = new MaterialEffects(fastMuHCAL);
0114 
0115   if (fastCalo.exists("ECALResponseScaling")) {
0116     ecalCorrection =
0117         std::make_unique<KKCorrectionFactors>(fastCalo.getParameter<edm::ParameterSet>("ECALResponseScaling"));
0118   }
0119 }
0120 
0121 void CalorimetryManager::clean() {
0122   EBMapping_.clear();
0123   EEMapping_.clear();
0124   HMapping_.clear();
0125   ESMapping_.clear();
0126   muonSimTracks.clear();
0127   savedMuonSimTracks.clear();
0128 }
0129 
0130 CalorimetryManager::~CalorimetryManager() {
0131   if (myCalorimeter_)
0132     delete myCalorimeter_;
0133   if (myHDResponse_)
0134     delete myHDResponse_;
0135 
0136   if (theMuonEcalEffects)
0137     delete theMuonEcalEffects;
0138   if (theMuonHcalEffects)
0139     delete theMuonHcalEffects;
0140 
0141   if (theProfile)
0142     delete theProfile;
0143 
0144   if (theHFShowerLibrary)
0145     delete theHFShowerLibrary;
0146 }
0147 
0148 void CalorimetryManager::reconstruct(RandomEngineAndDistribution const* random) {
0149   if (!evtsToDebug_.empty()) {
0150     std::vector<unsigned int>::const_iterator itcheck =
0151         find(evtsToDebug_.begin(), evtsToDebug_.end(), mySimEvent->id().event());
0152     debug_ = (itcheck != evtsToDebug_.end());
0153     if (debug_)
0154       mySimEvent->print();
0155   }
0156 
0157   initialize(random);
0158 
0159   LogInfo("FastCalorimetry") << "Reconstructing " << (int)mySimEvent->nTracks() << " tracks." << std::endl;
0160   for (int fsimi = 0; fsimi < (int)mySimEvent->nTracks(); ++fsimi) {
0161     FSimTrack& myTrack = mySimEvent->track(fsimi);
0162 
0163     reconstructTrack(myTrack, random);
0164   }  // particle loop
0165 
0166 }  // reconstruct
0167 
0168 void CalorimetryManager::initialize(RandomEngineAndDistribution const* random) {
0169   // Clear the content of the calorimeters
0170   if (!initialized_) {
0171     theHFShowerLibrary->SetRandom(random);
0172 
0173     // Check if the preshower is really available
0174     if (simulatePreshower_ && !myCalorimeter_->preshowerPresent()) {
0175       edm::LogWarning("CalorimetryManager")
0176           << " WARNING: The preshower simulation has been turned on; but no preshower geometry is available "
0177           << std::endl;
0178       edm::LogWarning("CalorimetryManager") << " Disabling the preshower simulation " << std::endl;
0179       simulatePreshower_ = false;
0180     }
0181 
0182     initialized_ = true;
0183   }
0184   clean();
0185 }
0186 
0187 void CalorimetryManager::reconstructTrack(FSimTrack& myTrack, RandomEngineAndDistribution const* random) {
0188   int pid = abs(myTrack.type());
0189 
0190   if (debug_) {
0191     LogInfo("FastCalorimetry") << " ===> pid = " << pid << std::endl;
0192   }
0193 
0194   // Check that the particle hasn't decayed
0195   if (myTrack.noEndVertex()) {
0196     // Simulate energy smearing for photon and electrons
0197     float charge_ = (float)(myTrack.charge());
0198     if (pid == 11 || pid == 22) {
0199       if (myTrack.onEcal())
0200         EMShowerSimulation(myTrack, random);
0201       else if (myTrack.onVFcal()) {
0202         if (useShowerLibrary) {
0203           theHFShowerLibrary->recoHFShowerLibrary(myTrack);
0204           myHDResponse_->correctHF(myTrack.hcalEntrance().e(), abs(myTrack.type()));
0205           updateHCAL(theHFShowerLibrary->getHitsMap(), myTrack.id());
0206         } else
0207           reconstructHCAL(myTrack, random);
0208       }
0209     }  // electron or photon
0210     else if (pid == 13 || pid == 1000024 || (pid > 1000100 && pid < 1999999 && fabs(charge_) > 0.001)) {
0211       MuonMipSimulation(myTrack, random);
0212     }
0213     // Simulate energy smearing for hadrons (i.e., everything
0214     // but muons... and SUSY particles that deserve a special
0215     // treatment.
0216     else if (pid < 1000000) {
0217       if (myTrack.onHcal() || myTrack.onVFcal()) {
0218         if (optionHDSim_ == 0)
0219           reconstructHCAL(myTrack, random);
0220         else
0221           HDShowerSimulation(myTrack, random);
0222       }
0223     }  // pid < 1000000
0224   }    // myTrack.noEndVertex()
0225 }
0226 
0227 // Simulation of electromagnetic showers in PS, ECAL, HCAL
0228 void CalorimetryManager::EMShowerSimulation(const FSimTrack& myTrack, RandomEngineAndDistribution const* random) {
0229   std::vector<const RawParticle*> thePart;
0230   double X0depth;
0231 
0232   if (debug_) {
0233     LogInfo("FastCalorimetry") << " EMShowerSimulation " << myTrack << std::endl;
0234   }
0235 
0236   // The Particle at ECAL entrance
0237   myPart = myTrack.ecalEntrance();
0238 
0239   // protection against infinite loop.
0240   if (myTrack.type() == 22 && myPart.e() < 0.055)
0241     return;
0242 
0243   // Barrel or Endcap ?
0244   int onEcal = myTrack.onEcal();
0245   int onHcal = myTrack.onHcal();
0246   int onLayer1 = myTrack.onLayer1();
0247   int onLayer2 = myTrack.onLayer2();
0248 
0249   // The entrance in ECAL
0250   XYZPoint ecalentrance = myPart.vertex().Vect();
0251 
0252   // The preshower
0253   PreshowerHitMaker* myPreshower = nullptr;
0254   if (simulatePreshower_ && (onLayer1 || onLayer2)) {
0255     XYZPoint layer1entrance, layer2entrance;
0256     XYZVector dir1, dir2;
0257     if (onLayer1) {
0258       layer1entrance = XYZPoint(myTrack.layer1Entrance().vertex().Vect());
0259       dir1 = XYZVector(myTrack.layer1Entrance().Vect().Unit());
0260     }
0261     if (onLayer2) {
0262       layer2entrance = XYZPoint(myTrack.layer2Entrance().vertex().Vect());
0263       dir2 = XYZVector(myTrack.layer2Entrance().Vect().Unit());
0264     }
0265     myPreshower =
0266         new PreshowerHitMaker(myCalorimeter_, layer1entrance, dir1, layer2entrance, dir2, aLandauGenerator, random);
0267     myPreshower->setMipEnergy(mipValues_[0], mipValues_[1]);
0268   }
0269 
0270   // The ECAL Properties
0271   EMECALShowerParametrization showerparam(myCalorimeter_->ecalProperties(onEcal),
0272                                           myCalorimeter_->hcalProperties(onHcal),
0273                                           myCalorimeter_->layer1Properties(onLayer1),
0274                                           myCalorimeter_->layer2Properties(onLayer2),
0275                                           theCoreIntervals_,
0276                                           theTailIntervals_,
0277                                           RCFactor_,
0278                                           RTFactor_);
0279 
0280   // Photons : create an e+e- pair
0281   if (myTrack.type() == 22) {
0282     // Depth for the first e+e- pair creation (in X0)
0283     X0depth = -log(random->flatShoot()) * (9. / 7.);
0284 
0285     // Initialization
0286     double eMass = 0.000510998902;
0287     double xe = 0;
0288     double xm = eMass / myPart.e();
0289     double weight = 0.;
0290 
0291     // Generate electron energy between emass and eGamma-emass
0292     do {
0293       xe = random->flatShoot() * (1. - 2. * xm) + xm;
0294       weight = 1. - 4. / 3. * xe * (1. - xe);
0295     } while (weight < random->flatShoot());
0296 
0297     // Protection agains infinite loop in Famos Shower
0298     if (myPart.e() * xe < 0.055 || myPart.e() * (1. - xe) < 0.055) {
0299       if (myPart.e() > 0.055)
0300         thePart.push_back(&myPart);
0301 
0302     } else {
0303       myElec = (myPart.momentum()) * xe;
0304       myPosi = (myPart.momentum()) * (1. - xe);
0305       myElec.setVertex(myPart.vertex());
0306       myPosi.setVertex(myPart.vertex());
0307       thePart.push_back(&myElec);
0308       thePart.push_back(&myPosi);
0309     }
0310     // Electrons
0311   } else {
0312     X0depth = 0.;
0313     if (myPart.e() > 0.055)
0314       thePart.push_back(&myPart);
0315   }
0316 
0317   // After the different protections, this shouldn't happen.
0318   if (thePart.empty()) {
0319     if (myPreshower == nullptr)
0320       return;
0321     delete myPreshower;
0322     return;
0323   }
0324 
0325   // find the most energetic particle
0326   double maxEnergy = -1.;
0327   for (unsigned ip = 0; ip < thePart.size(); ++ip)
0328     if (thePart[ip]->e() > maxEnergy)
0329       maxEnergy = thePart[ip]->e();
0330 
0331   // Initialize the Grid in ECAL
0332   int size = gridSize_;
0333   if (maxEnergy > 100)
0334     size = 11;
0335 
0336   EMShower theShower(random, aGammaGenerator, &showerparam, &thePart, nullptr, nullptr, bFixedLength_);
0337 
0338   double maxShower = theShower.getMaximumOfShower();
0339   if (maxShower > 20.)
0340     maxShower = 2.;  // simple pivot-searching protection
0341 
0342   double depth((X0depth + maxShower) * myCalorimeter_->ecalProperties(onEcal)->radLenIncm());
0343   XYZPoint meanShower = ecalentrance + myPart.Vect().Unit() * depth;
0344 
0345   // The closest crystal
0346   DetId pivot(myCalorimeter_->getClosestCell(meanShower, true, onEcal == 1));
0347 
0348   if (pivot.subdetId() == 0) {  // further protection against avbsence of pivot
0349     edm::LogWarning("CalorimetryManager")
0350         << "Pivot for egamma  e = " << myTrack.hcalEntrance().e() << " is not found at depth " << depth
0351         << " and meanShower coordinates = " << meanShower << std::endl;
0352     if (myPreshower)
0353       delete myPreshower;
0354     return;
0355   }
0356 
0357   EcalHitMaker myGrid(myCalorimeter_, ecalentrance, pivot, onEcal, size, 0, random);
0358   //                                             ^^^^
0359   //                                         for EM showers
0360   myGrid.setPulledPadSurvivalProbability(pulledPadSurvivalProbability_);
0361   myGrid.setCrackPadSurvivalProbability(crackPadSurvivalProbability_);
0362 
0363   //maximumdepth dependence of the radiusfactorbehindpreshower
0364   //First tuning: Shilpi Jain (Mar-Apr 2010); changed after tuning - Feb-July - Shilpi Jain
0365   /* **************
0366      myGrid.setRadiusFactor(radiusFactor_);
0367      if(onLayer1 || onLayer2)
0368      {
0369      float b               = radiusPreshowerCorrections_[0];
0370      float a               = radiusFactor_*( 1.+radiusPreshowerCorrections_[1]*radiusPreshowerCorrections_[0] );
0371      float maxdepth        = X0depth+theShower.getMaximumOfShower();
0372      float newRadiusFactor = radiusFactor_;
0373      if(myPart.e()<=250.)
0374      {
0375      newRadiusFactor = a/(1.+b*maxdepth); 
0376      }
0377      myGrid.setRadiusFactor(newRadiusFactor);
0378      }
0379      else // otherwise use the normal radius factor
0380      {
0381      myGrid.setRadiusFactor(radiusFactor_);
0382      }
0383      ************** */
0384   if (myTrack.onEcal() == 2)  // if on EE
0385   {
0386     if ((onLayer1 || onLayer2) && myPart.e() <= 250.) {
0387       double maxdepth = X0depth + theShower.getMaximumOfShower();
0388       double newRadiusFactor = radiusFactorEE_ * aTerm / (1. + bTerm * maxdepth);
0389       myGrid.setRadiusFactor(newRadiusFactor);
0390     } else  // otherwise use the normal radius factor
0391     {
0392       myGrid.setRadiusFactor(radiusFactorEE_);
0393     }
0394   }     //if(myTrack.onEcal() == 2)
0395   else  // else if on EB
0396   {
0397     myGrid.setRadiusFactor(radiusFactorEB_);
0398   }
0399   //(end of) changed after tuning - Feb-July - Shilpi Jain
0400 
0401   myGrid.setPreshowerPresent(simulatePreshower_);
0402 
0403   // The shower simulation
0404   myGrid.setTrackParameters(myPart.Vect().Unit(), X0depth, myTrack);
0405 
0406   if (myPreshower)
0407     theShower.setPreshower(myPreshower);
0408 
0409   HcalHitMaker myHcalHitMaker(myGrid, (unsigned)0);
0410 
0411   theShower.setGrid(&myGrid);
0412   theShower.setHcal(&myHcalHitMaker);
0413   theShower.compute();
0414 
0415   // calculate the total simulated energy for this particle
0416   float simE = 0;
0417   for (const auto& mapIterator : myGrid.getHits()) {
0418     simE += mapIterator.second;
0419   }
0420 
0421   auto scale = ecalCorrection
0422                    ? ecalCorrection->getScale(myTrack.ecalEntrance().e(), std::abs(myTrack.ecalEntrance().eta()), simE)
0423                    : 1.;
0424 
0425   // Save the hits !
0426   updateECAL(myGrid.getHits(), onEcal, myTrack.id(), scale);
0427 
0428   // Now fill the HCAL hits
0429   updateHCAL(myHcalHitMaker.getHits(), myTrack.id());
0430 
0431   // delete the preshower
0432   if (myPreshower != nullptr) {
0433     updatePreshower(myPreshower->getHits(), myTrack.id());
0434     delete myPreshower;
0435   }
0436 }
0437 
0438 void CalorimetryManager::reconstructHCAL(const FSimTrack& myTrack, RandomEngineAndDistribution const* random) {
0439   int hit;
0440   int pid = abs(myTrack.type());
0441   if (debug_) {
0442     LogInfo("FastCalorimetry") << " reconstructHCAL " << myTrack << std::endl;
0443   }
0444 
0445   XYZTLorentzVector trackPosition;
0446   if (myTrack.onHcal()) {
0447     trackPosition = myTrack.hcalEntrance().vertex();
0448     hit = myTrack.onHcal() - 1;
0449   } else {
0450     trackPosition = myTrack.vfcalEntrance().vertex();
0451     hit = 2;
0452   }
0453 
0454   double pathEta = trackPosition.eta();
0455   double pathPhi = trackPosition.phi();
0456 
0457   double EGen = myTrack.hcalEntrance().e();
0458   double emeas = 0.;
0459 
0460   float charge_ = (float)myTrack.charge();
0461   if (pid == 13 || pid == 1000024 || (pid > 1000100 && pid < 1999999 && fabs(charge_) > 0.001)) {
0462     emeas = myHDResponse_->responseHCAL(0, EGen, pathEta, 2, random);  // 2=muon
0463     if (debug_)
0464       LogInfo("FastCalorimetry") << "CalorimetryManager::reconstructHCAL - MUON !!!" << std::endl;
0465   } else if (pid == 22 || pid == 11) {
0466     emeas = myHDResponse_->responseHCAL(0, EGen, pathEta, 0, random);  // last par. = 0 = e/gamma
0467     if (debug_)
0468       LogInfo("FastCalorimetry") << "CalorimetryManager::reconstructHCAL - e/gamma !!!" << std::endl;
0469   } else {
0470     emeas = myHDResponse_->getHCALEnergyResponse(EGen, hit, random);
0471   }
0472 
0473   if (debug_)
0474     LogInfo("FastCalorimetry") << "CalorimetryManager::reconstructHCAL - on-calo "
0475                                << "  eta = " << pathEta << "  phi = " << pathPhi << "  Egen = " << EGen
0476                                << "  Emeas = " << emeas << std::endl;
0477 
0478   if (emeas > 0.) {
0479     DetId cell = myCalorimeter_->getClosestCell(trackPosition.Vect(), false, false);
0480     double tof =
0481         (((HcalGeometry*)(myCalorimeter_->getHcalGeometry()))->getPosition(cell).mag()) / 29.98;  //speed of light
0482     CaloHitID current_id(cell.rawId(), tof, myTrack.id());
0483     std::map<CaloHitID, float> hitMap;
0484     hitMap[current_id] = emeas;
0485     updateHCAL(hitMap, myTrack.id());
0486   }
0487 }
0488 
0489 void CalorimetryManager::HDShowerSimulation(const FSimTrack& myTrack, RandomEngineAndDistribution const* random) {  //,
0490   // const edm::ParameterSet& fastCalo){
0491 
0492   theHFShowerLibrary->SetRandom(random);
0493 
0494   //  TimeMe t(" FASTEnergyReconstructor::HDShower");
0495   const XYZTLorentzVector& moment = myTrack.momentum();
0496 
0497   if (debug_)
0498     LogInfo("FastCalorimetry") << "CalorimetryManager::HDShowerSimulation - track param." << std::endl
0499                                << "  eta = " << moment.eta() << std::endl
0500                                << "  phi = " << moment.phi() << std::endl
0501                                << "   et = " << moment.Et() << std::endl
0502                                << "   e  = " << myTrack.hcalEntrance().e() << std::endl;
0503 
0504   if (debug_) {
0505     LogInfo("FastCalorimetry") << " HDShowerSimulation " << myTrack << std::endl;
0506   }
0507 
0508   int hit;
0509 
0510   XYZTLorentzVector trackPosition;
0511   if (myTrack.onEcal()) {
0512     trackPosition = myTrack.ecalEntrance().vertex();
0513     hit = myTrack.onEcal() - 1;  //
0514     myPart = myTrack.ecalEntrance();
0515   } else if (myTrack.onVFcal()) {
0516     trackPosition = myTrack.vfcalEntrance().vertex();
0517     hit = 2;
0518     myPart = myTrack.vfcalEntrance();
0519   } else {
0520     LogInfo("FastCalorimetry") << " The particle is not in the acceptance " << std::endl;
0521     return;
0522   }
0523 
0524   // int onHCAL = hit + 1; - specially for myCalorimeter->hcalProperties(onHCAL)
0525   // (below) to get VFcal properties ...
0526   int onHCAL = hit + 1;
0527   int onECAL = myTrack.onEcal();
0528 
0529   double pathEta = trackPosition.eta();
0530   double pathPhi = trackPosition.phi();
0531 
0532   double eint = moment.e();
0533   double eGen = myTrack.hcalEntrance().e();
0534 
0535   double emeas = 0.;
0536   double pmip = myHDResponse_->getMIPfraction(eGen, pathEta);
0537 
0538   //===========================================================================
0539   if (eGen > 0.) {
0540     // ECAL and HCAL properties to get
0541     HDShowerParametrization theHDShowerparam(
0542         myCalorimeter_->ecalProperties(onECAL), myCalorimeter_->hcalProperties(onHCAL), myHSParameters_);
0543 
0544     //Making ECAL Grid (and segments calculation)
0545     XYZPoint caloentrance;
0546     XYZVector direction;
0547     if (myTrack.onEcal()) {
0548       caloentrance = myTrack.ecalEntrance().vertex().Vect();
0549       direction = myTrack.ecalEntrance().Vect().Unit();
0550     } else if (myTrack.onHcal()) {
0551       caloentrance = myTrack.hcalEntrance().vertex().Vect();
0552       direction = myTrack.hcalEntrance().Vect().Unit();
0553     } else {
0554       caloentrance = myTrack.vfcalEntrance().vertex().Vect();
0555       direction = myTrack.vfcalEntrance().Vect().Unit();
0556     }
0557 
0558     if (debug_)
0559       LogInfo("FastCalorimetry") << "CalorimetryManager::HDShowerSimulation - on-calo 1 " << std::endl
0560                                  << "  onEcal    = " << myTrack.onEcal() << std::endl
0561                                  << "  onHcal    = " << myTrack.onHcal() << std::endl
0562                                  << "  onVFcal   = " << myTrack.onVFcal() << std::endl
0563                                  << "  position  = " << caloentrance << std::endl;
0564 
0565     DetId pivot;
0566     if (myTrack.onEcal()) {
0567       pivot = myCalorimeter_->getClosestCell(caloentrance, true, myTrack.onEcal() == 1);
0568     } else if (myTrack.onHcal()) {
0569       pivot = myCalorimeter_->getClosestCell(caloentrance, false, false);
0570     }
0571 
0572     EcalHitMaker myGrid(
0573         myCalorimeter_, caloentrance, pivot, pivot.null() ? 0 : myTrack.onEcal(), hdGridSize_, 1, random);
0574     // 1=HAD shower
0575 
0576     myGrid.setTrackParameters(direction, 0, myTrack);
0577     // Build the FAMOS HCAL
0578     HcalHitMaker myHcalHitMaker(myGrid, (unsigned)1);
0579 
0580     // Shower simulation
0581     bool status = false;
0582     int mip = 2;
0583     // Use HFShower for HF
0584     if (!myTrack.onEcal() && !myTrack.onHcal()) {
0585       // Warning : We give here the particle energy with the response
0586       //           but without the resolution/gaussian smearing
0587       //           For HF, the resolution is due to the PE statistic
0588 
0589       if (useShowerLibrary) {
0590         theHFShowerLibrary->recoHFShowerLibrary(myTrack);
0591         status = true;
0592       } else {
0593         HFShower theShower(random, &theHDShowerparam, &myGrid, &myHcalHitMaker, onECAL, eGen);
0594         //           eGen);
0595         //           e); // PV Warning : temporarly set the energy to the generated E
0596 
0597         status = theShower.compute();
0598       }
0599     } else {
0600       if (hdSimMethod_ == 0) {
0601         HDShower theShower(random, &theHDShowerparam, &myGrid, &myHcalHitMaker, onECAL, eGen, pmip);
0602         status = theShower.compute();
0603         mip = theShower.getmip();
0604       } else if (hdSimMethod_ == 1) {
0605         HDRShower theShower(random, &theHDShowerparam, &myGrid, &myHcalHitMaker, onECAL, eGen);
0606         status = theShower.computeShower();
0607         mip = 2;
0608       } else if (hdSimMethod_ == 2) {
0609         //dynamically loading a corresponding profile by the particle type
0610         int particleType = myTrack.type();
0611         theProfile = thePiKProfile;
0612         if (particleType == -2212)
0613           theProfile = theAntiProtonProfile;
0614         else if (particleType == 2212)
0615           theProfile = theProtonProfile;
0616 
0617         //input variables for GflashHadronShowerProfile
0618         int showerType = 99 + myTrack.onEcal();
0619         double globalTime = 150.0;  // a temporary reference hit time in nanosecond
0620         float charge = (float)(myTrack.charge());
0621         Gflash3Vector gfpos(trackPosition.X(), trackPosition.Y(), trackPosition.Z());
0622         Gflash3Vector gfmom(moment.X(), moment.Y(), moment.Z());
0623 
0624         theProfile->initialize(showerType, eGen, globalTime, charge, gfpos, gfmom);
0625         theProfile->loadParameters();
0626         theProfile->hadronicParameterization();
0627 
0628         //make hits
0629         std::vector<GflashHit>& gflashHitList = theProfile->getGflashHitList();
0630         std::vector<GflashHit>::const_iterator spotIter = gflashHitList.begin();
0631         std::vector<GflashHit>::const_iterator spotIterEnd = gflashHitList.end();
0632 
0633         Gflash::CalorimeterNumber whichCalor = Gflash::kNULL;
0634 
0635         for (; spotIter != spotIterEnd; spotIter++) {
0636           double pathLength = theProfile->getGflashShowino()->getPathLengthAtShower() +
0637                               (30 * 100 / eGen) * (spotIter->getTime() - globalTime);
0638 
0639           double currentDepth = std::max(0.0, pathLength - theProfile->getGflashShowino()->getPathLengthOnEcal());
0640 
0641           //find the the showino position at the currentDepth
0642           GflashTrajectoryPoint trajectoryPoint;
0643           theProfile->getGflashShowino()->getHelix()->getGflashTrajectoryPoint(trajectoryPoint, pathLength);
0644           Gflash3Vector positionAtCurrentDepth = trajectoryPoint.getPosition();
0645           //find radial distrance
0646           Gflash3Vector lateralDisplacement = positionAtCurrentDepth - spotIter->getPosition() / CLHEP::cm;
0647           double rShower = lateralDisplacement.r();
0648           double azimuthalAngle = lateralDisplacement.phi();
0649 
0650           whichCalor = Gflash::getCalorimeterNumber(positionAtCurrentDepth);
0651 
0652           if (whichCalor == Gflash::kESPM || whichCalor == Gflash::kENCA) {
0653             bool statusPad = myGrid.getPads(currentDepth, true);
0654             if (!statusPad)
0655               continue;
0656             myGrid.setSpotEnergy(1.2 * spotIter->getEnergy() / CLHEP::GeV);
0657             myGrid.addHit(rShower / Gflash::intLength[Gflash::kESPM], azimuthalAngle, 0);
0658           } else if (whichCalor == Gflash::kHB || whichCalor == Gflash::kHE) {
0659             bool setHDdepth = myHcalHitMaker.setDepth(currentDepth, true);
0660             if (!setHDdepth)
0661               continue;
0662             myHcalHitMaker.setSpotEnergy(1.4 * spotIter->getEnergy() / CLHEP::GeV);
0663             myHcalHitMaker.addHit(rShower / Gflash::intLength[Gflash::kHB], azimuthalAngle, 0);
0664           }
0665         }
0666         status = true;
0667       } else {
0668         edm::LogInfo("FastSimulationCalorimetry") << " SimMethod " << hdSimMethod_ << " is NOT available ";
0669       }
0670     }
0671 
0672     if (status) {
0673       // Here to switch between simple formulae and parameterized response
0674       if (optionHDSim_ == 1) {
0675         emeas = myHDResponse_->getHCALEnergyResponse(eGen, hit, random);
0676       } else {                                                               // optionHDsim == 2
0677         emeas = myHDResponse_->responseHCAL(mip, eGen, pathEta, 1, random);  // 1=hadron
0678       }
0679 
0680       double correction = emeas / eGen;
0681 
0682       // RespCorrP factors (ECAL and HCAL separately) calculation
0683       respCorr(eint);
0684 
0685       if (debug_)
0686         LogInfo("FastCalorimetry") << "CalorimetryManager::HDShowerSimulation - on-calo 2" << std::endl
0687                                    << "   eta  = " << pathEta << std::endl
0688                                    << "   phi  = " << pathPhi << std::endl
0689                                    << "  Egen  = " << eGen << std::endl
0690                                    << " Emeas  = " << emeas << std::endl
0691                                    << "  corr  = " << correction << std::endl
0692                                    << "   mip  = " << mip << std::endl;
0693 
0694       if (myTrack.onEcal() > 0) {
0695         // Save ECAL hits
0696         updateECAL(myGrid.getHits(), onECAL, myTrack.id(), correction * ecorr);
0697       }
0698 
0699       // Save HCAL hits
0700       if (myTrack.onVFcal() && useShowerLibrary) {
0701         myHDResponse_->correctHF(eGen, abs(myTrack.type()));
0702         updateHCAL(theHFShowerLibrary->getHitsMap(), myTrack.id());
0703       } else
0704         updateHCAL(myHcalHitMaker.getHits(), myTrack.id(), correction * hcorr);
0705 
0706     } else {  // shower simulation failed
0707       if (myTrack.onHcal() || myTrack.onVFcal()) {
0708         DetId cell = myCalorimeter_->getClosestCell(trackPosition.Vect(), false, false);
0709         double tof =
0710             (((HcalGeometry*)(myCalorimeter_->getHcalGeometry()))->getPosition(cell).mag()) / 29.98;  //speed of light
0711         CaloHitID current_id(cell.rawId(), tof, myTrack.id());
0712         std::map<CaloHitID, float> hitMap;
0713         hitMap[current_id] = emeas;
0714         updateHCAL(hitMap, myTrack.id());
0715         if (debug_)
0716           LogInfo("FastCalorimetry") << " HCAL simple cell " << cell.rawId() << " added    E = " << emeas << std::endl;
0717       }
0718     }
0719 
0720   }  // e > 0. ...
0721 
0722   if (debug_)
0723     LogInfo("FastCalorimetry") << std::endl << " FASTEnergyReconstructor::HDShowerSimulation  finished " << std::endl;
0724 }
0725 
0726 void CalorimetryManager::MuonMipSimulation(const FSimTrack& myTrack, RandomEngineAndDistribution const* random) {
0727   //  TimeMe t(" FASTEnergyReconstructor::HDShower");
0728   XYZTLorentzVector moment = myTrack.momentum();
0729 
0730   // Backward compatibility behaviour
0731   if (!theMuonHcalEffects) {
0732     savedMuonSimTracks.push_back(myTrack);
0733 
0734     if (myTrack.onHcal() || myTrack.onVFcal())
0735       reconstructHCAL(myTrack, random);
0736 
0737     return;
0738   }
0739 
0740   if (debug_)
0741     LogInfo("FastCalorimetry") << "CalorimetryManager::MuonMipSimulation - track param." << std::endl
0742                                << "  eta = " << moment.eta() << std::endl
0743                                << "  phi = " << moment.phi() << std::endl
0744                                << "   et = " << moment.Et() << std::endl;
0745 
0746   XYZTLorentzVector trackPosition;
0747   if (myTrack.onEcal()) {
0748     trackPosition = myTrack.ecalEntrance().vertex();
0749     myPart = myTrack.ecalEntrance();
0750   } else if (myTrack.onVFcal()) {
0751     trackPosition = myTrack.vfcalEntrance().vertex();
0752     myPart = myTrack.vfcalEntrance();
0753   } else {
0754     LogInfo("FastCalorimetry") << " The particle is not in the acceptance " << std::endl;
0755     return;
0756   }
0757 
0758   // int onHCAL = hit + 1; - specially for myCalorimeter->hcalProperties(onHCAL)
0759   // (below) to get VFcal properties ...
0760   // not needed ?
0761   //  int onHCAL = hit + 1;
0762   int onECAL = myTrack.onEcal();
0763 
0764   //===========================================================================
0765 
0766   // ECAL and HCAL properties to get
0767 
0768   //Making ECAL Grid (and segments calculation)
0769   XYZPoint caloentrance;
0770   XYZVector direction;
0771   if (myTrack.onEcal()) {
0772     caloentrance = myTrack.ecalEntrance().vertex().Vect();
0773     direction = myTrack.ecalEntrance().Vect().Unit();
0774   } else if (myTrack.onHcal()) {
0775     caloentrance = myTrack.hcalEntrance().vertex().Vect();
0776     direction = myTrack.hcalEntrance().Vect().Unit();
0777   } else {
0778     caloentrance = myTrack.vfcalEntrance().vertex().Vect();
0779     direction = myTrack.vfcalEntrance().Vect().Unit();
0780   }
0781 
0782   DetId pivot;
0783   if (myTrack.onEcal()) {
0784     pivot = myCalorimeter_->getClosestCell(caloentrance, true, myTrack.onEcal() == 1);
0785   } else if (myTrack.onHcal()) {
0786     pivot = myCalorimeter_->getClosestCell(caloentrance, false, false);
0787   }
0788 
0789   EcalHitMaker myGrid(myCalorimeter_, caloentrance, pivot, pivot.null() ? 0 : myTrack.onEcal(), hdGridSize_, 0, random);
0790   // 0 =EM shower -> Unit = X0
0791 
0792   myGrid.setTrackParameters(direction, 0, myTrack);
0793 
0794   // Now get the path in the Preshower, ECAL and HCAL along a straight line extrapolation
0795   // but only those in the ECAL are used
0796 
0797   const std::vector<CaloSegment>& segments = myGrid.getSegments();
0798   unsigned nsegments = segments.size();
0799 
0800   int ifirstHcal = -1;
0801   int ilastEcal = -1;
0802 
0803   EnergyLossSimulator* energyLossECAL = (theMuonEcalEffects) ? theMuonEcalEffects->energyLossSimulator() : nullptr;
0804   //  // Muon brem in ECAL
0805   //  MuonBremsstrahlungSimulator* muonBremECAL = 0;
0806   //  if (theMuonEcalEffects) muonBremECAL = theMuonEcalEffects->muonBremsstrahlungSimulator();
0807 
0808   for (unsigned iseg = 0; iseg < nsegments && ifirstHcal < 0; ++iseg) {
0809     // in the ECAL, there are two types of segments: PbWO4 and GAP
0810     float segmentSizeinX0 = segments[iseg].X0length();
0811 
0812     // Martijn - insert your computations here
0813     float energy = 0.0;
0814     if (segmentSizeinX0 > 0.001 && segments[iseg].material() == CaloSegment::PbWO4) {
0815       // The energy loss simulator
0816       float charge = (float)(myTrack.charge());
0817       RawParticle p = rawparticle::makeMuon(charge < 0, moment, trackPosition);
0818       ParticlePropagator theMuon(p, nullptr, nullptr, mySimEvent->theTable());
0819       if (energyLossECAL) {
0820         energyLossECAL->updateState(theMuon, segmentSizeinX0, random);
0821         energy = energyLossECAL->deltaMom().E();
0822         moment -= energyLossECAL->deltaMom();
0823       }
0824     }
0825     // that's all for ECAL, Florian
0826     // Save the hit only if it is a crystal
0827     if (segments[iseg].material() == CaloSegment::PbWO4) {
0828       myGrid.getPads(segments[iseg].sX0Entrance() + segmentSizeinX0 * 0.5);
0829       myGrid.setSpotEnergy(energy);
0830       myGrid.addHit(0., 0.);
0831       ilastEcal = iseg;
0832     }
0833     // Check for end of loop:
0834     if (segments[iseg].material() == CaloSegment::HCAL) {
0835       ifirstHcal = iseg;
0836     }
0837   }
0838 
0839   // Build the FAMOS HCAL
0840   HcalHitMaker myHcalHitMaker(myGrid, (unsigned)2);
0841   // float mipenergy=0.1;
0842   // Create the helix with the stepping helix propagator
0843   // to add a hit, just do
0844   // myHcalHitMaker.setSpotEnergy(mipenergy);
0845   // math::XYZVector hcalEntrance;
0846   // if(ifirstHcal>=0) hcalEntrance=segments[ifirstHcal].entrance();
0847   // myHcalHitMaker.addHit(hcalEntrance);
0848   ///
0849   /////
0850   ////// TEMPORARY First attempt to include HCAL (with straight-line extrapolation):
0851   int ilastHcal = -1;
0852   float mipenergy = 0.0;
0853 
0854   EnergyLossSimulator* energyLossHCAL = (theMuonHcalEffects) ? theMuonHcalEffects->energyLossSimulator() : nullptr;
0855   //  // Muon Brem effect
0856   //  MuonBremsstrahlungSimulator* muonBremHCAL = 0;
0857   //  if (theMuonHcalEffects) muonBremHCAL = theMuonHcalEffects->muonBremsstrahlungSimulator();
0858 
0859   if (ifirstHcal > 0 && energyLossHCAL) {
0860     for (unsigned iseg = ifirstHcal; iseg < nsegments; ++iseg) {
0861       float segmentSizeinX0 = segments[iseg].X0length();
0862       if (segments[iseg].material() == CaloSegment::HCAL) {
0863         ilastHcal = iseg;
0864         if (segmentSizeinX0 > 0.001) {
0865           // The energy loss simulator
0866           float charge = (float)(myTrack.charge());
0867           RawParticle p = rawparticle::makeMuon(charge < 0, moment, trackPosition);
0868           ParticlePropagator theMuon(p, nullptr, nullptr, mySimEvent->theTable());
0869           energyLossHCAL->updateState(theMuon, segmentSizeinX0, random);
0870           mipenergy = energyLossHCAL->deltaMom().E();
0871           moment -= energyLossHCAL->deltaMom();
0872           myHcalHitMaker.setSpotEnergy(mipenergy);
0873           myHcalHitMaker.addHit(segments[iseg].entrance());
0874         }
0875       }
0876     }
0877   }
0878 
0879   // Copy the muon SimTrack (Only for Energy loss)
0880   FSimTrack muonTrack(myTrack);
0881   if (energyLossHCAL && ilastHcal >= 0) {
0882     math::XYZVector hcalExit = segments[ilastHcal].exit();
0883     muonTrack.setTkPosition(hcalExit);
0884     muonTrack.setTkMomentum(moment);
0885   } else if (energyLossECAL && ilastEcal >= 0) {
0886     math::XYZVector ecalExit = segments[ilastEcal].exit();
0887     muonTrack.setTkPosition(ecalExit);
0888     muonTrack.setTkMomentum(moment);
0889   }  // else just leave tracker surface position and momentum...
0890 
0891   muonSimTracks.push_back(muonTrack);
0892 
0893   // no need to change below this line
0894   std::map<CaloHitID, float>::const_iterator mapitr;
0895   std::map<CaloHitID, float>::const_iterator endmapitr;
0896   if (myTrack.onEcal() > 0) {
0897     // Save ECAL hits
0898     updateECAL(myGrid.getHits(), onECAL, myTrack.id());
0899   }
0900 
0901   // Save HCAL hits
0902   updateHCAL(myHcalHitMaker.getHits(), myTrack.id());
0903 
0904   if (debug_)
0905     LogInfo("FastCalorimetry") << std::endl << " FASTEnergyReconstructor::MipShowerSimulation  finished " << std::endl;
0906 }
0907 
0908 void CalorimetryManager::readParameters(const edm::ParameterSet& fastCalo) {
0909   edm::ParameterSet ECALparameters = fastCalo.getParameter<edm::ParameterSet>("ECAL");
0910 
0911   evtsToDebug_ = fastCalo.getUntrackedParameter<std::vector<unsigned int> >("EvtsToDebug", std::vector<unsigned>());
0912   debug_ = fastCalo.getUntrackedParameter<bool>("Debug");
0913 
0914   bFixedLength_ = ECALparameters.getParameter<bool>("bFixedLength");
0915 
0916   gridSize_ = ECALparameters.getParameter<int>("GridSize");
0917   spotFraction_ = ECALparameters.getParameter<double>("SpotFraction");
0918   pulledPadSurvivalProbability_ = ECALparameters.getParameter<double>("FrontLeakageProbability");
0919   crackPadSurvivalProbability_ = ECALparameters.getParameter<double>("GapLossProbability");
0920   theCoreIntervals_ = ECALparameters.getParameter<std::vector<double> >("CoreIntervals");
0921   theTailIntervals_ = ECALparameters.getParameter<std::vector<double> >("TailIntervals");
0922 
0923   RCFactor_ = ECALparameters.getParameter<double>("RCFactor");
0924   RTFactor_ = ECALparameters.getParameter<double>("RTFactor");
0925   //changed after tuning - Feb-July - Shilpi Jain
0926   //  radiusFactor_ = ECALparameters.getParameter<double>("RadiusFactor");
0927   radiusFactorEE_ = ECALparameters.getParameter<double>("RadiusFactorEE");
0928   radiusFactorEB_ = ECALparameters.getParameter<double>("RadiusFactorEB");
0929   //(end of) changed after tuning - Feb-July - Shilpi Jain
0930   radiusPreshowerCorrections_ = ECALparameters.getParameter<std::vector<double> >("RadiusPreshowerCorrections");
0931   aTerm = 1. + radiusPreshowerCorrections_[1] * radiusPreshowerCorrections_[0];
0932   bTerm = radiusPreshowerCorrections_[0];
0933   mipValues_ = ECALparameters.getParameter<std::vector<double> >("MipsinGeV");
0934   simulatePreshower_ = ECALparameters.getParameter<bool>("SimulatePreshower");
0935 
0936   if (gridSize_ < 1)
0937     gridSize_ = 7;
0938   if (pulledPadSurvivalProbability_ < 0. || pulledPadSurvivalProbability_ > 1)
0939     pulledPadSurvivalProbability_ = 1.;
0940   if (crackPadSurvivalProbability_ < 0. || crackPadSurvivalProbability_ > 1)
0941     crackPadSurvivalProbability_ = 0.9;
0942 
0943   LogInfo("FastCalorimetry") << " Fast ECAL simulation parameters " << std::endl;
0944   LogInfo("FastCalorimetry") << " =============================== " << std::endl;
0945   if (simulatePreshower_)
0946     LogInfo("FastCalorimetry") << " The preshower is present " << std::endl;
0947   else
0948     LogInfo("FastCalorimetry") << " The preshower is NOT present " << std::endl;
0949   LogInfo("FastCalorimetry") << " Grid Size : " << gridSize_ << std::endl;
0950   if (spotFraction_ > 0.)
0951     LogInfo("FastCalorimetry") << " Spot Fraction : " << spotFraction_ << std::endl;
0952   else {
0953     LogInfo("FastCalorimetry") << " Core of the shower " << std::endl;
0954     for (unsigned ir = 0; ir < theCoreIntervals_.size() / 2; ++ir) {
0955       LogInfo("FastCalorimetry") << " r < " << theCoreIntervals_[ir * 2] << " R_M : " << theCoreIntervals_[ir * 2 + 1]
0956                                  << "        ";
0957     }
0958     LogInfo("FastCalorimetry") << std::endl;
0959 
0960     LogInfo("FastCalorimetry") << " Tail of the shower " << std::endl;
0961     for (unsigned ir = 0; ir < theTailIntervals_.size() / 2; ++ir) {
0962       LogInfo("FastCalorimetry") << " r < " << theTailIntervals_[ir * 2] << " R_M : " << theTailIntervals_[ir * 2 + 1]
0963                                  << "        ";
0964     }
0965     //changed after tuning - Feb-July - Shilpi Jain
0966     LogInfo("FastCalorimetry") << "Radius correction factors:  EB & EE " << radiusFactorEB_ << " : " << radiusFactorEE_
0967                                << std::endl;
0968     //(end of) changed after tuning - Feb-July - Shilpi Jain
0969     LogInfo("FastCalorimetry") << std::endl;
0970     if (mipValues_.size() > 2) {
0971       LogInfo("FastCalorimetry") << "Improper number of parameters for the preshower ; using 95keV" << std::endl;
0972       mipValues_.clear();
0973       mipValues_.resize(2, 0.000095);
0974     }
0975   }
0976 
0977   LogInfo("FastCalorimetry") << " FrontLeakageProbability : " << pulledPadSurvivalProbability_ << std::endl;
0978   LogInfo("FastCalorimetry") << " GapLossProbability : " << crackPadSurvivalProbability_ << std::endl;
0979 
0980   // RespCorrP: p (momentum), ECAL and HCAL corrections = f(p)
0981   edm::ParameterSet CalorimeterParam = fastCalo.getParameter<edm::ParameterSet>("CalorimeterProperties");
0982 
0983   rsp = CalorimeterParam.getParameter<std::vector<double> >("RespCorrP");
0984   LogInfo("FastCalorimetry") << " RespCorrP (rsp) size " << rsp.size() << std::endl;
0985 
0986   if (rsp.size() % 3 != 0) {
0987     LogInfo("FastCalorimetry") << " RespCorrP size is wrong -> no corrections applied !!!" << std::endl;
0988 
0989     p_knots.push_back(14000.);
0990     k_e.push_back(1.);
0991     k_h.push_back(1.);
0992   } else {
0993     for (unsigned i = 0; i < rsp.size(); i += 3) {
0994       LogInfo("FastCalorimetry") << "i = " << i / 3 << "   p = " << rsp[i] << "   k_e(p) = " << rsp[i + 1]
0995                                  << "   k_e(p) = " << rsp[i + 2] << std::endl;
0996 
0997       p_knots.push_back(rsp[i]);
0998       k_e.push_back(rsp[i + 1]);
0999       k_h.push_back(rsp[i + 2]);
1000     }
1001   }
1002 
1003   //FR
1004   edm::ParameterSet HCALparameters = fastCalo.getParameter<edm::ParameterSet>("HCAL");
1005   optionHDSim_ = HCALparameters.getParameter<int>("SimOption");
1006   hdGridSize_ = HCALparameters.getParameter<int>("GridSize");
1007   hdSimMethod_ = HCALparameters.getParameter<int>("SimMethod");
1008   //RF
1009 
1010   EcalDigitizer_ = ECALparameters.getUntrackedParameter<bool>("Digitizer", false);
1011   HcalDigitizer_ = HCALparameters.getUntrackedParameter<bool>("Digitizer", false);
1012   samplingHBHE_ = HCALparameters.getParameter<std::vector<double> >("samplingHBHE");
1013   samplingHF_ = HCALparameters.getParameter<std::vector<double> >("samplingHF");
1014   samplingHO_ = HCALparameters.getParameter<std::vector<double> >("samplingHO");
1015   ietaShiftHB_ = HCALparameters.getParameter<int>("ietaShiftHB");
1016   ietaShiftHE_ = HCALparameters.getParameter<int>("ietaShiftHE");
1017   ietaShiftHF_ = HCALparameters.getParameter<int>("ietaShiftHF");
1018   ietaShiftHO_ = HCALparameters.getParameter<int>("ietaShiftHO");
1019   timeShiftHB_ = HCALparameters.getParameter<std::vector<double> >("timeShiftHB");
1020   timeShiftHE_ = HCALparameters.getParameter<std::vector<double> >("timeShiftHE");
1021   timeShiftHF_ = HCALparameters.getParameter<std::vector<double> >("timeShiftHF");
1022   timeShiftHO_ = HCALparameters.getParameter<std::vector<double> >("timeShiftHO");
1023 
1024   // FastHFShowerLibrary
1025   edm::ParameterSet m_HS = fastCalo.getParameter<edm::ParameterSet>("HFShowerLibrary");
1026   useShowerLibrary = m_HS.getUntrackedParameter<bool>("useShowerLibrary", false);
1027   useCorrectionSL = m_HS.getUntrackedParameter<bool>("useCorrectionSL", false);
1028 }
1029 
1030 void CalorimetryManager::respCorr(double p) {
1031   int sizeP = p_knots.size();
1032 
1033   if (sizeP <= 1) {
1034     ecorr = 1.;
1035     hcorr = 1.;
1036   } else {
1037     int ip = -1;
1038     for (int i = 0; i < sizeP; i++) {
1039       if (p < p_knots[i]) {
1040         ip = i;
1041         break;
1042       }
1043     }
1044     if (ip == 0) {
1045       ecorr = k_e[0];
1046       hcorr = k_h[0];
1047     } else {
1048       if (ip == -1) {
1049         ecorr = k_e[sizeP - 1];
1050         hcorr = k_h[sizeP - 1];
1051       } else {
1052         double x1 = p_knots[ip - 1];
1053         double x2 = p_knots[ip];
1054         double y1 = k_e[ip - 1];
1055         double y2 = k_e[ip];
1056 
1057         ecorr = (y1 + (y2 - y1) * (p - x1) / (x2 - x1));
1058 
1059         y1 = k_h[ip - 1];
1060         y2 = k_h[ip];
1061         hcorr = (y1 + (y2 - y1) * (p - x1) / (x2 - x1));
1062       }
1063     }
1064   }
1065 
1066   if (debug_)
1067     LogInfo("FastCalorimetry") << " p, ecorr, hcorr = " << p << " " << ecorr << "  " << hcorr << std::endl;
1068 }
1069 
1070 void CalorimetryManager::updateECAL(const std::map<CaloHitID, float>& hitMap, int onEcal, int trackID, float corr) {
1071   std::map<CaloHitID, float>::const_iterator mapitr;
1072   std::map<CaloHitID, float>::const_iterator endmapitr = hitMap.end();
1073   if (onEcal == 1) {
1074     EBMapping_.reserve(EBMapping_.size() + hitMap.size());
1075     endmapitr = hitMap.end();
1076     for (mapitr = hitMap.begin(); mapitr != endmapitr; ++mapitr) {
1077       //correct energy
1078       float energy = mapitr->second;
1079       energy *= corr;
1080 
1081       //make finalized CaloHitID
1082       CaloHitID current_id(mapitr->first.unitID(), mapitr->first.timeSlice(), trackID);
1083 
1084       EBMapping_.push_back(std::pair<CaloHitID, float>(current_id, energy));
1085     }
1086   } else if (onEcal == 2) {
1087     EEMapping_.reserve(EEMapping_.size() + hitMap.size());
1088     endmapitr = hitMap.end();
1089     for (mapitr = hitMap.begin(); mapitr != endmapitr; ++mapitr) {
1090       //correct energy
1091       float energy = mapitr->second;
1092       energy *= corr;
1093 
1094       //make finalized CaloHitID
1095       CaloHitID current_id(mapitr->first.unitID(), mapitr->first.timeSlice(), trackID);
1096 
1097       EEMapping_.push_back(std::pair<CaloHitID, float>(current_id, energy));
1098     }
1099   }
1100 }
1101 
1102 void CalorimetryManager::updateHCAL(const std::map<CaloHitID, float>& hitMap, int trackID, float corr) {
1103   std::vector<double> hfcorrEm = myHDResponse_->getCorrHFem();
1104   std::vector<double> hfcorrHad = myHDResponse_->getCorrHFhad();
1105   std::map<CaloHitID, float>::const_iterator mapitr;
1106   std::map<CaloHitID, float>::const_iterator endmapitr = hitMap.end();
1107   HMapping_.reserve(HMapping_.size() + hitMap.size());
1108   for (mapitr = hitMap.begin(); mapitr != endmapitr; ++mapitr) {
1109     //correct energy
1110     float energy = mapitr->second;
1111     energy *= corr;
1112 
1113     float time = mapitr->first.timeSlice();
1114     //put energy into uncalibrated state for digitizer && correct timing
1115     if (HcalDigitizer_) {
1116       HcalDetId hdetid = HcalDetId(mapitr->first.unitID());
1117       if (hdetid.subdetId() == HcalBarrel) {
1118         energy /= samplingHBHE_[hdetid.ietaAbs() - 1];  //re-convert to GeV
1119         time = timeShiftHB_[hdetid.ietaAbs() - ietaShiftHB_];
1120       } else if (hdetid.subdetId() == HcalEndcap) {
1121         energy /= samplingHBHE_[hdetid.ietaAbs() - 1];  //re-convert to GeV
1122         time = timeShiftHE_[hdetid.ietaAbs() - ietaShiftHE_];
1123       } else if (hdetid.subdetId() == HcalForward) {
1124         if (useShowerLibrary) {
1125           if (useCorrectionSL) {
1126             if (hdetid.depth() == 1 or hdetid.depth() == 3)
1127               energy *= hfcorrEm[hdetid.ietaAbs() - ietaShiftHF_];
1128             if (hdetid.depth() == 2 or hdetid.depth() == 4)
1129               energy *= hfcorrHad[hdetid.ietaAbs() - ietaShiftHF_];
1130           }
1131         } else {
1132           if (hdetid.depth() == 1 or hdetid.depth() == 3)
1133             energy *= samplingHF_[0];
1134           if (hdetid.depth() == 2 or hdetid.depth() == 4)
1135             energy *= samplingHF_[1];
1136           time = timeShiftHF_[hdetid.ietaAbs() - ietaShiftHF_];
1137         }
1138       } else if (hdetid.subdetId() == HcalOuter) {
1139         energy /= samplingHO_[hdetid.ietaAbs() - 1];
1140         time = timeShiftHO_[hdetid.ietaAbs() - ietaShiftHO_];
1141       }
1142     }
1143 
1144     //make finalized CaloHitID
1145     CaloHitID current_id(mapitr->first.unitID(), time, trackID);
1146     HMapping_.push_back(std::pair<CaloHitID, float>(current_id, energy));
1147   }
1148 }
1149 
1150 void CalorimetryManager::updatePreshower(const std::map<CaloHitID, float>& hitMap, int trackID, float corr) {
1151   std::map<CaloHitID, float>::const_iterator mapitr;
1152   std::map<CaloHitID, float>::const_iterator endmapitr = hitMap.end();
1153   ESMapping_.reserve(ESMapping_.size() + hitMap.size());
1154   for (mapitr = hitMap.begin(); mapitr != endmapitr; ++mapitr) {
1155     //correct energy
1156     float energy = mapitr->second;
1157     energy *= corr;
1158 
1159     //make finalized CaloHitID
1160     CaloHitID current_id(mapitr->first.unitID(), mapitr->first.timeSlice(), trackID);
1161 
1162     ESMapping_.push_back(std::pair<CaloHitID, float>(current_id, energy));
1163   }
1164 }
1165 
1166 void CalorimetryManager::loadFromEcalBarrel(edm::PCaloHitContainer& c) const {
1167   c.reserve(c.size() + EBMapping_.size());
1168   for (unsigned i = 0; i < EBMapping_.size(); i++) {
1169     c.push_back(PCaloHit(EBDetId::unhashIndex(EBMapping_[i].first.unitID()),
1170                          EBMapping_[i].second,
1171                          EBMapping_[i].first.timeSlice(),
1172                          EBMapping_[i].first.trackID()));
1173   }
1174 }
1175 
1176 void CalorimetryManager::loadFromEcalEndcap(edm::PCaloHitContainer& c) const {
1177   c.reserve(c.size() + EEMapping_.size());
1178   for (unsigned i = 0; i < EEMapping_.size(); i++) {
1179     c.push_back(PCaloHit(EEDetId::unhashIndex(EEMapping_[i].first.unitID()),
1180                          EEMapping_[i].second,
1181                          EEMapping_[i].first.timeSlice(),
1182                          EEMapping_[i].first.trackID()));
1183   }
1184 }
1185 
1186 void CalorimetryManager::loadFromHcal(edm::PCaloHitContainer& c) const {
1187   c.reserve(c.size() + HMapping_.size());
1188   for (unsigned i = 0; i < HMapping_.size(); i++) {
1189     c.push_back(PCaloHit(DetId(HMapping_[i].first.unitID()),
1190                          HMapping_[i].second,
1191                          HMapping_[i].first.timeSlice(),
1192                          HMapping_[i].first.trackID()));
1193   }
1194 }
1195 
1196 void CalorimetryManager::loadFromPreshower(edm::PCaloHitContainer& c) const {
1197   c.reserve(c.size() + ESMapping_.size());
1198   for (unsigned i = 0; i < ESMapping_.size(); i++) {
1199     c.push_back(PCaloHit(ESMapping_[i].first.unitID(),
1200                          ESMapping_[i].second,
1201                          ESMapping_[i].first.timeSlice(),
1202                          ESMapping_[i].first.trackID()));
1203   }
1204 }
1205 
1206 // The main danger in this method is to screw up to relationships between particles
1207 // So, the muon FSimTracks created by FSimEvent.cc are simply to be updated
1208 void CalorimetryManager::loadMuonSimTracks(edm::SimTrackContainer& muons) const {
1209   unsigned size = muons.size();
1210   for (unsigned i = 0; i < size; ++i) {
1211     int id = muons[i].trackId();
1212     if (!(abs(muons[i].type()) == 13 || abs(muons[i].type()) == 1000024 ||
1213           (abs(muons[i].type()) > 1000100 && abs(muons[i].type()) < 1999999)))
1214       continue;
1215     // identify the corresponding muon in the local collection
1216 
1217     std::vector<FSimTrack>::const_iterator itcheck =
1218         find_if(muonSimTracks.begin(), muonSimTracks.end(), FSimTrackEqual(id));
1219     if (itcheck != muonSimTracks.end()) {
1220       muons[i].setTkPosition(itcheck->trackerSurfacePosition());
1221       muons[i].setTkMomentum(itcheck->trackerSurfaceMomentum());
1222     }
1223   }
1224 }
1225 
1226 void CalorimetryManager::harvestMuonSimTracks(edm::SimTrackContainer& c) const {
1227   c.reserve(int(0.2 * muonSimTracks.size() + 0.2 * savedMuonSimTracks.size() + 0.5));
1228   for (const auto& track : muonSimTracks) {
1229     if (track.momentum().perp2() > 1.0 && fabs(track.momentum().eta()) < 3.0 && track.isGlobal())
1230       c.push_back(track);
1231   }
1232   for (const auto& track : savedMuonSimTracks) {
1233     if (track.momentum().perp2() > 1.0 && fabs(track.momentum().eta()) < 3.0 && track.isGlobal())
1234       c.push_back(track);
1235   }
1236   c.shrink_to_fit();
1237 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!