Backport of "Fix for FastSim decays of exotic-descendent SM particles decaying outside pipe" to 10_2_X

Configuration/ProcessModifiers/python/fastSimFixLongLivedBug_cff.py

@@ -0,0 +1,4 @@
+import FWCore.ParameterSet.Config as cms
+
+# Desiged to disable the bug in Run II samples that duplicates hits for long lived particles
+fastSimFixLongLivedBug =  cms.Modifier()

FastSimulation/SimplifiedGeometryPropagator/interface/Decayer.h

@@ -50,10 +50,11 @@ namespace fastsim
             \param engine The Random Engine.
         */
         void decay(const Particle & particle, std::vector<std::unique_ptr<Particle> > & secondaries, CLHEP::HepRandomEngine & engine) const;
-
+        void setfixLongLivedBug(bool disable){ fixLongLivedBug_ = disable; };            
         private:    
         std::unique_ptr<Pythia8::Pythia> pythia_;  //!< Instance of pythia
         std::unique_ptr<gen::P8RndmEngine> pythiaRandomEngine_;  //!< Instance of pythia Random Engine
+        bool fixLongLivedBug_;
     };
 }
 #endif

FastSimulation/SimplifiedGeometryPropagator/interface/ParticleManager.h

@@ -56,7 +56,8 @@ namespace fastsim {
             double deltaRchargedMother,
             const ParticleFilter & particleFilter,
             std::vector<SimTrack> & simTracks,
-            std::vector<SimVertex> & simVertices);
+            std::vector<SimVertex> & simVertices,
+            bool fixLongLivedBug);
 
         //! Default destructor.
         ~ParticleManager();
@@ -142,15 +143,18 @@ namespace fastsim {
         double lengthUnitConversionFactor2_;  //!< Convert pythia unis to cm^2 (FastSim standard)
         double timeUnitConversionFactor_;  //!< Convert pythia unis to ns (FastSim standard)
         std::vector<std::unique_ptr<Particle> > particleBuffer_;  //!< The vector of all secondaries that are not yet propagated in the event.
+        bool fixLongLivedBug_;
     };
 }
 
-inline bool isExotic(int pdgid_) {
+inline bool isExotic(bool fixLongLivedBug, int pdgid_) {
   unsigned int pdgid = std::abs(pdgid_);
   return ((pdgid >= 1000000 && pdgid < 4000000 && pdgid != 3000022) ||  // SUSY, R-hadron, and technicolor particles
           pdgid == 17 ||                                                // 4th generation lepton
           pdgid == 34 ||                                                // W-prime
-          pdgid == 37);                                                 // charged Higgs
+          pdgid == 37 ||                                                // charged Higgs  
+          (pdgid == 39 && fixLongLivedBug));                           // bulk graviton
+
 }
 
 #endif

FastSimulation/SimplifiedGeometryPropagator/plugins/FastSimProducer.cc

@@ -93,6 +93,8 @@ class FastSimProducer : public edm::stream::EDProducer<> {
     std::vector<std::unique_ptr<fastsim::InteractionModel> > interactionModels_;  //!< All defined interaction models
     std::map<std::string, fastsim::InteractionModel *> interactionModelMap_;  //!< Each interaction model has a unique name
     static const std::string MESSAGECATEGORY;  //!< Category of debugging messages ("FastSimulation")
+
+    bool fixLongLivedBug_;
 };
 
 const std::string FastSimProducer::MESSAGECATEGORY = "FastSimulation";
@@ -107,8 +109,10 @@ FastSimProducer::FastSimProducer(const edm::ParameterSet& iConfig)
     , _randomEngine(nullptr)
     , simulateCalorimetry(iConfig.getParameter<bool>("simulateCalorimetry"))
     , simulateMuons(iConfig.getParameter<bool>("simulateMuons"))
+    , fixLongLivedBug_(iConfig.getParameter<bool>("fixLongLivedBug"))
 {
-
+    // Fix the decayer for the long lived stuff
+    decayer_.setfixLongLivedBug(fixLongLivedBug_);
     //----------------
     // define interaction models
     //---------------
@@ -195,7 +199,8 @@ FastSimProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
                                             ,deltaRchargedMother_
                                             ,particleFilter_
                                             ,*simTracks_
-                                            ,*simVertices_);
+                                            ,*simVertices_
+                                            ,fixLongLivedBug_);
 
     //  Initialize the calorimeter geometry
     if(simulateCalorimetry)

FastSimulation/SimplifiedGeometryPropagator/python/fastSimProducer_cff.py

@@ -77,4 +77,9 @@
     MaterialEffectsForMuonsInECAL = MaterialEffectsForMuonsInECALBlock.MaterialEffectsForMuonsInECAL,
     MaterialEffectsForMuonsInHCAL = MaterialEffectsForMuonsInHCALBlock.MaterialEffectsForMuonsInHCAL,
     GFlash = FamosCalorimetryBlock.GFlash,
+    fixLongLivedBug = cms.bool(False),
 )
+
+from Configuration.ProcessModifiers.fastSimFixLongLivedBug_cff import fastSimFixLongLivedBug
+
+fastSimFixLongLivedBug.toModify(fastSimProducer, fixLongLivedBug = cms.bool(True))

FastSimulation/SimplifiedGeometryPropagator/src/Decayer.cc

@@ -18,7 +18,7 @@ fastsim::Decayer::Decayer()
     pythia_->settings.flag("PartonLevel:FSRinResonances",false);
     pythia_->settings.flag("ProcessLevel:resonanceDecays",false);
     pythia_->init();
-
+    fixLongLivedBug_ = false;
     // forbid all decays
     // (decays are allowed selectively in the decay function)
     Pythia8::ParticleData & pdt = pythia_->particleData;
@@ -39,9 +39,8 @@ fastsim::Decayer::decay(const Particle & particle,std::vector<std::unique_ptr<fa
     // inspired by method Pythia8Hadronizer::residualDecay() in GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc
     int pid = particle.pdgId();  
     // snip decay products of exotic particles or their children. These decay products are preserved from the event record.
-    // limitation: if exotic incurs heavy energy loss during propagation, the saved decay products could be too hard.
-
-    if (isExotic(pid) || isExotic(particle.getMotherPdgId())) {
+    // limitation: if exotic incurs heavy energy loss during propagation, the saved decay products could be too hard.    
+    if (isExotic(fixLongLivedBug_, pid) || isExotic(fixLongLivedBug_, particle.getMotherPdgId())) {
         return;
     }
 

FastSimulation/SimplifiedGeometryPropagator/src/ParticleManager.cc

@@ -21,7 +21,8 @@ fastsim::ParticleManager::ParticleManager(
     double deltaRchargedMother,
     const fastsim::ParticleFilter & particleFilter,
     std::vector<SimTrack> & simTracks,
-    std::vector<SimVertex> & simVertices)
+    std::vector<SimVertex> & simVertices,
+    bool fixLongLivedBug)
     : genEvent_(&genEvent)
     , genParticleIterator_(genEvent_->particles_begin())
     , genParticleEnd_(genEvent_->particles_end())
@@ -44,7 +45,7 @@ fastsim::ParticleManager::ParticleManager(
     , lengthUnitConversionFactor_(conversion_factor(genEvent_->length_unit(),HepMC::Units::LengthUnit::CM))
     , lengthUnitConversionFactor2_(lengthUnitConversionFactor_*lengthUnitConversionFactor_)
     , timeUnitConversionFactor_(lengthUnitConversionFactor_/fastsim::Constants::speedOfLight)
-
+    , fixLongLivedBug_(fixLongLivedBug)
 {
 
     // add the main vertex from the signal event to the simvertex collection
@@ -240,7 +241,6 @@ std::unique_ptr<fastsim::Particle> fastsim::ParticleManager::nextGenParticle()
         const HepMC::GenParticle & particle = **genParticleIterator_;
         const HepMC::GenVertex * productionVertex = particle.production_vertex();
         const HepMC::GenVertex * endVertex = particle.end_vertex();
-
         // skip incoming particles
         if(!productionVertex){
             continue;
@@ -251,22 +251,27 @@ std::unique_ptr<fastsim::Particle> fastsim::ParticleManager::nextGenParticle()
 
         // particles which do not descend from exotics must be produced within the beampipe
         int exoticRelativeId = 0;
-        if (productionVertex->position().perp2() * lengthUnitConversionFactor2_ > beamPipeRadius2_)  //
+        const bool producedWithinBeamPipe = productionVertex->position().perp2() * lengthUnitConversionFactor2_ < beamPipeRadius2_;
+        if (!producedWithinBeamPipe)  //
         {
             exoticRelativesChecker(productionVertex, exoticRelativeId, 0);
-            if (!isExotic(exoticRelativeId)) {
+            if (!isExotic(fixLongLivedBug_, exoticRelativeId)) {
                 continue;
             }
-        }
-
-
-        // particle must not decay before it reaches the beam pipe
-        if(endVertex && endVertex->position().perp2()*lengthUnitConversionFactor2_ < beamPipeRadius2_)
+        }	
+        const bool decayedWithinBeamPipe = endVertex && endVertex->position().perp2() * lengthUnitConversionFactor2_ < beamPipeRadius2_;
+        if(decayedWithinBeamPipe)
         {
             continue;
         }
 
-
+        // SM particles that descend from exotics and cross the beam pipe radius should make hits but not be decayed,
+        // by default it will duplicate FastSim hits for long lived particles 
+        // and so anything produced without activating fixLongLivedBug_ is physically wrong
+        if (fixLongLivedBug_ && producedWithinBeamPipe && !decayedWithinBeamPipe){
+          exoticRelativesChecker(productionVertex, exoticRelativeId, 0);
+        } 
+
         // make the particle
         std::unique_ptr<Particle> newParticle(
             new Particle(particle.pdg_id(),
@@ -279,7 +284,7 @@ std::unique_ptr<fastsim::Particle> fastsim::ParticleManager::nextGenParticle()
                              particle.momentum().z()*momentumUnitConversionFactor_,
                              particle.momentum().e()*momentumUnitConversionFactor_)));
         newParticle->setGenParticleIndex(genParticleIndex_);
-        if (isExotic(exoticRelativeId)) {
+        if (isExotic(fixLongLivedBug_, exoticRelativeId)) {
             newParticle->setMotherPdgId(exoticRelativeId);
         }
         // try to get the life time of the particle from the genEvent
@@ -313,14 +318,14 @@ std::unique_ptr<fastsim::Particle> fastsim::ParticleManager::nextGenParticle()
 void fastsim::ParticleManager::exoticRelativesChecker(const HepMC::GenVertex* originVertex,
                                                       int& exoticRelativeId_,
                                                       int ngendepth = 0) {
-  if (ngendepth > 99 || exoticRelativeId_ == -1 || isExotic(std::abs(exoticRelativeId_)))
+  if (ngendepth > 99 || exoticRelativeId_ == -1 || isExotic(fixLongLivedBug_, std::abs(exoticRelativeId_)))
     return;
   ngendepth += 1;
   std::vector<HepMC::GenParticle*>::const_iterator relativesIterator_ = originVertex->particles_in_const_begin();
   std::vector<HepMC::GenParticle*>::const_iterator relativesIteratorEnd_ = originVertex->particles_in_const_end();
   for (; relativesIterator_ != relativesIteratorEnd_; ++relativesIterator_) {
     const HepMC::GenParticle& genRelative = **relativesIterator_;
-    if (isExotic(std::abs(genRelative.pdg_id()))) {
+    if (isExotic(fixLongLivedBug_, std::abs(genRelative.pdg_id()))) {
       exoticRelativeId_ = genRelative.pdg_id();
       if (ngendepth == 100)
         exoticRelativeId_ = -1;