Merge pull request #36352 from cericeci/backportFastSim_10_6_X

Backport of "Fix for FastSim decays of exotic-descendent SM particles decaying outside pipe" to 10_6_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();
@@ -143,15 +144,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
     //---------------
@@ -194,7 +198,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;
@@ -40,7 +40,7 @@ fastsim::Decayer::decay(const Particle & particle,std::vector<std::unique_ptr<fa
     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())) {
+    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;
@@ -252,20 +252,26 @@ 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_;
+        // FastSim will not make hits out of particles that decay before reaching the beam pipe
+        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(),
@@ -278,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
@@ -312,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;