diff --git a/Configuration/StandardSequences/python/Digi_cff.py b/Configuration/StandardSequences/python/Digi_cff.py index 64f516b0b2a0d..69cec91601cef 100644 --- a/Configuration/StandardSequences/python/Digi_cff.py +++ b/Configuration/StandardSequences/python/Digi_cff.py @@ -28,11 +28,12 @@ # add updating the GEN information by default from Configuration.StandardSequences.Generator_cff import * from GeneratorInterface.Core.generatorSmeared_cfi import * +from SimGeneral.PileupInformation.genPUProtons_cfi import * doAllDigi = cms.Sequence(generatorSmeared*calDigi+muonDigi) -pdigi = cms.Sequence(generatorSmeared*fixGenInfo*cms.SequencePlaceholder("randomEngineStateProducer")*cms.SequencePlaceholder("mix")*doAllDigi*addPileupInfo) +pdigi = cms.Sequence(generatorSmeared*fixGenInfo*cms.SequencePlaceholder("randomEngineStateProducer")*cms.SequencePlaceholder("mix")*doAllDigi*addPileupInfo*genPUProtons) pdigi_valid = cms.Sequence(pdigi) -pdigi_nogen=cms.Sequence(generatorSmeared*cms.SequencePlaceholder("randomEngineStateProducer")*cms.SequencePlaceholder("mix")*doAllDigi*addPileupInfo) +pdigi_nogen=cms.Sequence(generatorSmeared*cms.SequencePlaceholder("randomEngineStateProducer")*cms.SequencePlaceholder("mix")*doAllDigi*addPileupInfo*genPUProtons) pdigi_valid_nogen=cms.Sequence(pdigi_nogen) from GeneratorInterface.HiGenCommon.HeavyIon_cff import * diff --git a/PhysicsTools/PatAlgos/python/slimming/MicroEventContent_cff.py b/PhysicsTools/PatAlgos/python/slimming/MicroEventContent_cff.py index 6e4ba1449cfe4..920a6025f1829 100644 --- a/PhysicsTools/PatAlgos/python/slimming/MicroEventContent_cff.py +++ b/PhysicsTools/PatAlgos/python/slimming/MicroEventContent_cff.py @@ -87,6 +87,7 @@ 'keep GenLumiInfoHeader_generator_*_*', 'keep GenLumiInfoProduct_*_*_*', 'keep GenEventInfoProduct_generator_*_*', + 'keep recoGenParticles_genPUProtons_*_*', 'keep *_slimmedGenJetsFlavourInfos_*_*', 'keep *_slimmedGenJets__*', 'keep *_slimmedGenJetsAK8__*', diff --git a/SimGeneral/Configuration/python/SimGeneral_EventContent_cff.py b/SimGeneral/Configuration/python/SimGeneral_EventContent_cff.py index 4be3b39d522c1..3001602b0113a 100644 --- a/SimGeneral/Configuration/python/SimGeneral_EventContent_cff.py +++ b/SimGeneral/Configuration/python/SimGeneral_EventContent_cff.py @@ -12,17 +12,20 @@ outputCommands = cms.untracked.vstring('keep CrossingFramePlaybackInfoNew_*_*_*', 'keep PileupSummaryInfos_*_*_*', 'keep int6stdbitsetstdpairs_*_AffectedAPVList_*', - 'keep int_*_bunchSpacing_*') + 'keep int_*_bunchSpacing_*', + 'keep *_genPUProtons_*_*') ) #RECO content SimGeneralRECO = cms.PSet( outputCommands = cms.untracked.vstring('keep PileupSummaryInfos_*_*_*', - 'keep int_*_bunchSpacing_*') + 'keep int_*_bunchSpacing_*', + 'keep *_genPUProtons_*_*') ) #AOD content SimGeneralAOD = cms.PSet( outputCommands = cms.untracked.vstring('keep PileupSummaryInfos_*_*_*', - 'keep int_*_bunchSpacing_*') + 'keep int_*_bunchSpacing_*', + 'keep *_genPUProtons_*_*') ) # mods for HGCAL diff --git a/SimGeneral/DataMixingModule/plugins/BuildFile.xml b/SimGeneral/DataMixingModule/plugins/BuildFile.xml index 17c874906783b..3ed775f2108fa 100644 --- a/SimGeneral/DataMixingModule/plugins/BuildFile.xml +++ b/SimGeneral/DataMixingModule/plugins/BuildFile.xml @@ -17,6 +17,7 @@ + diff --git a/SimGeneral/DataMixingModule/plugins/DataMixingModule.cc b/SimGeneral/DataMixingModule/plugins/DataMixingModule.cc index e211b32c5bb10..95d2c3de73832 100644 --- a/SimGeneral/DataMixingModule/plugins/DataMixingModule.cc +++ b/SimGeneral/DataMixingModule/plugins/DataMixingModule.cc @@ -24,6 +24,7 @@ // // #include "SimDataFormats/CrossingFrame/interface/CrossingFramePlaybackInfoNew.h" +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" #include "DataMixingModule.h" #include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h" @@ -235,6 +236,12 @@ namespace edm produces< int >("bunchSpacing"); produces(); + std::vector GenPUProtonsInputTags; + GenPUProtonsInputTags = ps.getParameter >("GenPUProtonsInputTags"); + for(std::vector::const_iterator it_InputTag = GenPUProtonsInputTags.begin(); + it_InputTag != GenPUProtonsInputTags.end(); ++it_InputTag) + produces< std::vector >( it_InputTag->label() ); + PUWorker_ = new DataMixingPileupCopy(ps, consumesCollector()); } diff --git a/SimGeneral/DataMixingModule/plugins/DataMixingPileupCopy.cc b/SimGeneral/DataMixingModule/plugins/DataMixingPileupCopy.cc index dd0fd74c6a755..e2e8daab768ad 100644 --- a/SimGeneral/DataMixingModule/plugins/DataMixingPileupCopy.cc +++ b/SimGeneral/DataMixingModule/plugins/DataMixingPileupCopy.cc @@ -12,6 +12,8 @@ #include "DataFormats/Common/interface/Handle.h" #include "DataFormats/Provenance/interface/Provenance.h" #include "DataFormats/Provenance/interface/BranchDescription.h" + +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" // // #include "DataMixingPileupCopy.h" @@ -38,6 +40,7 @@ namespace edm BunchSpacingInputTag_ = ps.getParameter("BunchSpacingInputTag"); CFPlaybackInputTag_ = ps.getParameter("CFPlaybackInputTag"); + GenPUProtonsInputTags_ = ps.getParameter >("GenPUProtonsInputTags"); // apparently, we don't need consumes from Secondary input stream //iC.consumes>(PileupInfoInputTag_); @@ -77,6 +80,17 @@ namespace edm bsStorage_=10000; } + // Gen. PU protons + std::shared_ptr > const> GenPUProtonsPTR; + for(std::vector::const_iterator it_InputTag = GenPUProtonsInputTags_.begin(); + it_InputTag != GenPUProtonsInputTags_.end(); ++it_InputTag){ + GenPUProtonsPTR = getProductByTag >( *ep, *it_InputTag , mcc); + if( GenPUProtonsPTR != nullptr ){ + GenPUProtons_.push_back( *( GenPUProtonsPTR->product() ) ); + GenPUProtons_labels_.push_back( it_InputTag->label() ); + } else edm::LogWarning("DataMixingPileupCopy") << "Missing product with label: " << ( *it_InputTag ).label(); + } + // Playback std::shared_ptr const> PlaybackPTR = getProductByTag(*ep,CFPlaybackInputTag_, mcc); @@ -105,7 +119,19 @@ namespace edm e.put(std::move(PSIVector)); e.put(std::move(bsInt),"bunchSpacing"); + // Gen. PU protons + for(size_t idx = 0; idx < GenPUProtons_.size(); ++idx){ + std::unique_ptr > GenPUProtons_ptr( new std::vector() ); + std::vector::const_iterator it_GenParticle = GenPUProtons_.at(idx).begin(); + std::vector::const_iterator it_GenParticle_end = GenPUProtons_.at(idx).end(); + for(; it_GenParticle != it_GenParticle_end; ++ it_GenParticle) GenPUProtons_ptr->push_back( *it_GenParticle ); + + e.put( std::move(GenPUProtons_ptr), GenPUProtons_labels_.at(idx) ); + } + // clear local storage after this event PileupSummaryStorage_.clear(); + GenPUProtons_.clear(); + GenPUProtons_labels_.clear(); } } //edm diff --git a/SimGeneral/DataMixingModule/plugins/DataMixingPileupCopy.h b/SimGeneral/DataMixingModule/plugins/DataMixingPileupCopy.h index 3eba5a75c8488..15729b07751e9 100644 --- a/SimGeneral/DataMixingModule/plugins/DataMixingPileupCopy.h +++ b/SimGeneral/DataMixingModule/plugins/DataMixingPileupCopy.h @@ -29,6 +29,9 @@ #include #include +namespace reco{ + class GenParticle; +} namespace edm { @@ -61,12 +64,16 @@ namespace edm edm::InputTag BunchSpacingInputTag_ ; // InputTag for bunch spacing int edm::InputTag CFPlaybackInputTag_ ; // InputTag for CrossingFrame Playback information + std::vector GenPUProtonsInputTags_ ; CrossingFramePlaybackInfoNew CrossingFramePlaybackStorage_; std::vector PileupSummaryStorage_; int bsStorage_; + std::vector GenPUProtons_labels_; + std::vector > GenPUProtons_; + // unsigned int eventId_; //=0 for signal, from 1-n for pileup events std::string label_; diff --git a/SimGeneral/DataMixingModule/python/mixOne_sim_on_sim_cfi.py b/SimGeneral/DataMixingModule/python/mixOne_sim_on_sim_cfi.py index 3f5d5505027f8..8aae886968617 100644 --- a/SimGeneral/DataMixingModule/python/mixOne_sim_on_sim_cfi.py +++ b/SimGeneral/DataMixingModule/python/mixOne_sim_on_sim_cfi.py @@ -36,6 +36,7 @@ PileupInfoInputTag = cms.InputTag("addPileupInfo"), BunchSpacingInputTag = cms.InputTag("addPileupInfo","bunchSpacing"), CFPlaybackInputTag = cms.InputTag("mix"), + GenPUProtonsInputTags = cms.VInputTag("genPUProtons"), # SistripLabelSig = cms.InputTag("simSiStripDigis","ZeroSuppressed"), # diff --git a/SimGeneral/DataMixingModule/python/mixOne_simraw_on_sim_cfi.py b/SimGeneral/DataMixingModule/python/mixOne_simraw_on_sim_cfi.py index 38c073cba1353..b9716f4728098 100644 --- a/SimGeneral/DataMixingModule/python/mixOne_simraw_on_sim_cfi.py +++ b/SimGeneral/DataMixingModule/python/mixOne_simraw_on_sim_cfi.py @@ -122,6 +122,7 @@ PileupInfoInputTag = cms.InputTag("addPileupInfo"), BunchSpacingInputTag = cms.InputTag("addPileupInfo","bunchSpacing"), CFPlaybackInputTag = cms.InputTag("mix"), + GenPUProtonsInputTags = cms.VInputTag("genPUProtons"), # SistripLabelSig = cms.InputTag("simSiStripDigis","ZeroSuppressed"), # diff --git a/SimGeneral/MixingModule/python/aliases_PreMix_cfi.py b/SimGeneral/MixingModule/python/aliases_PreMix_cfi.py index e95c34bbfcde4..9117512db8f8c 100644 --- a/SimGeneral/MixingModule/python/aliases_PreMix_cfi.py +++ b/SimGeneral/MixingModule/python/aliases_PreMix_cfi.py @@ -40,6 +40,12 @@ # ) #) +genPUProtons = cms.EDAlias( + mixData = cms.VPSet( + cms.PSet( type = cms.string('recoGenParticles') ) + ) +) + # no castor,pixel,strip digis in fastsim from Configuration.Eras.Modifier_fastSim_cff import fastSim fastSim.toModify(simCastorDigis, mix = None) diff --git a/SimGeneral/MixingModule/python/mixObjects_cfi.py b/SimGeneral/MixingModule/python/mixObjects_cfi.py index 8c98b63758a5f..df042027c7c49 100644 --- a/SimGeneral/MixingModule/python/mixObjects_cfi.py +++ b/SimGeneral/MixingModule/python/mixObjects_cfi.py @@ -119,7 +119,7 @@ mixSimVertices.input = cms.VInputTag(cms.InputTag("famosSimHits")) mixHepMCProducts = cms.PSet( - makeCrossingFrame = cms.untracked.bool(False), + makeCrossingFrame = cms.untracked.bool(True), input = cms.VInputTag(cms.InputTag("generatorSmeared"),cms.InputTag("generator")), type = cms.string('HepMCProduct') ) diff --git a/SimGeneral/PileupInformation/plugins/BuildFile.xml b/SimGeneral/PileupInformation/plugins/BuildFile.xml index 4f6a220ecdc96..85cd2568d2410 100644 --- a/SimGeneral/PileupInformation/plugins/BuildFile.xml +++ b/SimGeneral/PileupInformation/plugins/BuildFile.xml @@ -11,6 +11,7 @@ + diff --git a/SimGeneral/PileupInformation/plugins/GenPUProtonProducer.cc b/SimGeneral/PileupInformation/plugins/GenPUProtonProducer.cc new file mode 100755 index 0000000000000..d399098173618 --- /dev/null +++ b/SimGeneral/PileupInformation/plugins/GenPUProtonProducer.cc @@ -0,0 +1,255 @@ +/* \class GenPUProtonProducer + * + * Modification of GenParticleProducer. + * Saves final state protons from HepMC events in Crossing Frame, in the generator-particle format. + * + * Note: Use the option USER_CXXFLAGS=-DEDM_ML_DEBUG with SCRAM in order to enable debug messages. + * + * March 9, 2017 : Initial version. + * March 14, 2017 : Updated debug messages. + * July 27, 2017 : Removed extra loop initially inherited from GenParticleProducer. + * August 17, 2017 : Replaced std::auto_ptr with std::unique_ptr. + * September 6, 2017 : Updated module to edm::global::EDProducer with ConvertParticle as RunCache following GenParticleProducer. + * + */ + +#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" +#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h" +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" +#include "FWCore/Utilities/interface/EDMException.h" + +#include +#include +#include + +namespace { + + class ConvertParticle { + public: + static constexpr int PDGCacheMax = 32768; + static constexpr double mmToCm = 0.1; + + ConvertParticle() : + abortOnUnknownPDGCode_( true ), + initialized_( false ), + chargeP_( PDGCacheMax, 0 ), chargeM_( PDGCacheMax, 0 ) {}; + + ConvertParticle(bool abortOnUnknownPDGCode) : + abortOnUnknownPDGCode_( abortOnUnknownPDGCode ), + initialized_( false ), + chargeP_( PDGCacheMax, 0 ), chargeM_( PDGCacheMax, 0 ) {}; + + ~ConvertParticle() {}; + + bool initialized() const { return initialized_; } + + void init(HepPDT::ParticleDataTable const& pdt) { + if( !initialized_ ) { + for( HepPDT::ParticleDataTable::const_iterator p = pdt.begin(); p != pdt.end(); ++p ) { + HepPDT::ParticleID const& id = p->first; + int pdgId = id.pid(), apdgId = std::abs( pdgId ); + int q3 = id.threeCharge(); + if ( apdgId < PDGCacheMax && pdgId > 0 ) { + chargeP_[ apdgId ] = q3; + chargeM_[ apdgId ] = -q3; + } else if ( apdgId < PDGCacheMax ) { + chargeP_[ apdgId ] = -q3; + chargeM_[ apdgId ] = q3; + } else { + chargeMap_.emplace( pdgId, q3); + chargeMap_.emplace( -pdgId, -q3); + } + } + initialized_ = true; + } + } + + bool operator() (reco::GenParticle& cand, HepMC::GenParticle const* part) const { + reco::Candidate::LorentzVector p4( part->momentum() ); + int pdgId = part->pdg_id(); + cand.setThreeCharge( chargeTimesThree( pdgId ) ); + cand.setPdgId( pdgId ); + cand.setStatus( part->status() ); + cand.setP4( p4 ); + cand.setCollisionId(0); + HepMC::GenVertex const* v = part->production_vertex(); + if ( v != nullptr ) { + HepMC::ThreeVector vtx = v->point3d(); + reco::Candidate::Point vertex( vtx.x() * mmToCm, vtx.y() * mmToCm, vtx.z() * mmToCm ); + cand.setVertex( vertex ); + } else { + cand.setVertex( reco::Candidate::Point( 0, 0, 0 ) ); + } + return true; + } + + private: + bool abortOnUnknownPDGCode_; + bool initialized_; + std::vector chargeP_, chargeM_; + std::unordered_map chargeMap_; + + int chargeTimesThree( int id ) const { + if( std::abs( id ) < PDGCacheMax ) + return id > 0 ? chargeP_[ id ] : chargeM_[ - id ]; + + auto f = chargeMap_.find( id ); + if ( f == chargeMap_.end() ) { + if ( abortOnUnknownPDGCode_ ) + throw edm::Exception( edm::errors::LogicError ) << "invalid PDG id: " << id << std::endl; + else + return HepPDT::ParticleID(id).threeCharge(); + } + return f->second; + } + + }; + + class SelectProton { + public: + bool operator() ( HepMC::GenParticle const* part, double minPz ) const { + bool selection = ( ( !part->end_vertex() && part->status() == 1 ) && + ( part->pdg_id() == 2212 ) && ( TMath::Abs( part->momentum().pz() ) >= minPz ) ); + return selection; + } + }; + +} // Anonymous namespace + +#include "FWCore/Framework/interface/global/EDProducer.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "DataFormats/Candidate/interface/CandidateFwd.h" +#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h" + +#include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h" +#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" +#include "PhysicsTools/HepMCCandAlgos/interface/MCTruthHelper.h" +#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h" +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" + +namespace edm { class ParameterSet; } + +class GenPUProtonProducer : public edm::global::EDProducer > { + public: + GenPUProtonProducer( const edm::ParameterSet & ); + ~GenPUProtonProducer() override; + + void produce( edm::StreamID, edm::Event& e, const edm::EventSetup&) const override; + std::shared_ptr globalBeginRun(const edm::Run&, const edm::EventSetup&) const override; + void globalEndRun(edm::Run const&, edm::EventSetup const&) const override {}; + + private: + edm::EDGetTokenT > mixToken_; + + bool abortOnUnknownPDGCode_; + std::vector bunchList_; + double minPz_; + SelectProton select_; + +}; + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "DataFormats/Common/interface/Handle.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/Run.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Utilities/interface/EDMException.h" + +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" +#include "SimDataFormats/CrossingFrame/interface/MixCollection.h" +#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h" + +#include + +using namespace edm; +using namespace reco; +using namespace std; +using namespace HepMC; + +GenPUProtonProducer::GenPUProtonProducer( const ParameterSet & cfg ) : + abortOnUnknownPDGCode_( cfg.getUntrackedParameter( "abortOnUnknownPDGCode", true ) ), + bunchList_( cfg.getParameter >( "bunchCrossingList" ) ), + minPz_( cfg.getParameter( "minPz" ) ) +{ + produces(); + mixToken_ = consumes >( InputTag(cfg.getParameter( "mix" ),"generatorSmeared") ); +} + +GenPUProtonProducer::~GenPUProtonProducer() { } + +std::shared_ptr GenPUProtonProducer::globalBeginRun(const Run&, const EventSetup& es) const { + ESHandle pdt; + es.getData( pdt ); + auto convert_ptr = std::make_shared( abortOnUnknownPDGCode_ ); + if( !convert_ptr->initialized() ) convert_ptr->init( *pdt ); + + return convert_ptr; +} + +void GenPUProtonProducer::produce( StreamID, Event& evt, const EventSetup& es ) const { + + size_t totalSize = 0; + size_t npiles = 1; + + Handle > cf; + evt.getByToken(mixToken_,cf); + std::unique_ptr > cfhepmcprod( new MixCollection( cf.product() ) ); + npiles = cfhepmcprod->size(); + + LogDebug("GenPUProtonProducer") << " Number of pile-up events : " << npiles << endl; + + for(size_t icf = 0; icf < npiles; ++icf){ + LogDebug("GenPUProtonProducer") << "CF " << icf << " size : " << cfhepmcprod->getObject(icf).GetEvent()->particles_size() << endl; + totalSize += cfhepmcprod->getObject(icf).GetEvent()->particles_size(); + } + LogDebug("GenPUProtonProducer") << "Total size : " << totalSize << endl; + + // Initialise containers + auto candsPtr = std::make_unique(); + GenParticleCollection& cands = *candsPtr; + + // Loop over pile-up events + ConvertParticle const& convertParticle_ = *runCache( evt.getRun().index() ); + + MixCollection::MixItr mixHepMC_itr; + unsigned int total_number_of_protons = 0; + size_t idx_mix = 0; + // Fill collection + for( mixHepMC_itr = cfhepmcprod->begin() ; mixHepMC_itr != cfhepmcprod->end() ; ++mixHepMC_itr, ++idx_mix ){ + int bunch = mixHepMC_itr.bunch(); + + if( find( bunchList_.begin(), bunchList_.end(), bunch ) != bunchList_.end() ){ + + auto event = (*mixHepMC_itr).GetEvent(); + + size_t num_particles = event->particles_size(); + + // Fill output collection + unsigned int number_of_protons = 0; + for( auto p = event->particles_begin() ; p != event->particles_end() ; ++p ) { + HepMC::GenParticle const* part = *p; + if( select_(part, minPz_) ) { + reco::GenParticle cand; + convertParticle_( cand, part ); + ++number_of_protons; + cands.push_back( cand ); + } + } + LogDebug("GenPUProtonProducer") << "Idx : " << idx_mix << " Bunch : " << bunch + << " Number of particles : " << num_particles + << " Number of protons : " << number_of_protons << endl; + + total_number_of_protons += number_of_protons; + } + } + LogDebug("GenPUProtonProducer") << "Total number of protons : " << total_number_of_protons << endl; + LogDebug("GenPUProtonProducer") << "Output collection size : " << cands.size() << endl; + + evt.put( std::move(candsPtr) ); +} + +#include "FWCore/Framework/interface/MakerMacros.h" + +DEFINE_FWK_MODULE( GenPUProtonProducer ); + diff --git a/SimGeneral/PileupInformation/python/genPUProtons_cfi.py b/SimGeneral/PileupInformation/python/genPUProtons_cfi.py new file mode 100644 index 0000000000000..875845cbab6b0 --- /dev/null +++ b/SimGeneral/PileupInformation/python/genPUProtons_cfi.py @@ -0,0 +1,7 @@ +import FWCore.ParameterSet.Config as cms + +genPUProtons = cms.EDProducer("GenPUProtonProducer", + mix = cms.string("mix"), + bunchCrossingList = cms.vint32(0), + minPz = cms.double( 2400. ) + )