| 1 |
|
| 2 |
#include <map>
|
| 3 |
#include <memory>
|
| 4 |
#include <string>
|
| 5 |
#include <utility>
|
| 6 |
#include <vector>
|
| 7 |
|
| 8 |
#include "DataFormats/Common/interface/Handle.h"
|
| 9 |
|
| 10 |
#include "FWCore/Framework/interface/EDProducer.h"
|
| 11 |
#include "FWCore/Framework/interface/Frameworkfwd.h"
|
| 12 |
#include "FWCore/Framework/interface/MakerMacros.h"
|
| 13 |
#include "FWCore/ParameterSet/interface/ParameterSet.h"
|
| 14 |
|
| 15 |
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
|
| 16 |
#include "SimDataFormats/Track/interface/SimTrack.h"
|
| 17 |
#include "SimDataFormats/Track/interface/SimTrackContainer.h"
|
| 18 |
#include "SimDataFormats/Vertex/interface/SimVertex.h"
|
| 19 |
#include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
|
| 20 |
|
| 21 |
//! Producer to reduce the number of charge in simulation
|
| 22 |
class HiEventTuner : public edm::EDProducer
|
| 23 |
{
|
| 24 |
public:
|
| 25 |
|
| 26 |
//! Contructor
|
| 27 |
explicit HiEventTuner(const edm::ParameterSet&);
|
| 28 |
|
| 29 |
//! Destructor
|
| 30 |
~HiEventTuner(){};
|
| 31 |
|
| 32 |
private:
|
| 33 |
|
| 34 |
static std::string MessageCategory_;
|
| 35 |
|
| 36 |
virtual void produce(edm::Event&, const edm::EventSetup&);
|
| 37 |
|
| 38 |
bool useMultipleHepMCLabels_;
|
| 39 |
std::vector<std::string> hepMCLabels_;
|
| 40 |
std::vector<edm::Handle<edm::HepMCProduct> > hepMCProducts_;
|
| 41 |
|
| 42 |
std::string simHitLabel_;
|
| 43 |
|
| 44 |
edm::Handle<edm::SimTrackContainer> simTracks_;
|
| 45 |
edm::Handle<edm::SimVertexContainer> simVertexes_;
|
| 46 |
|
| 47 |
std::auto_ptr<edm::SimTrackContainer> tunedSimTracks_;
|
| 48 |
std::auto_ptr<edm::SimVertexContainer> tunedSimVertexes_;
|
| 49 |
|
| 50 |
typedef std::multimap<std::size_t, std::size_t> IndexToIndexes;
|
| 51 |
|
| 52 |
IndexToIndexes trackIndexToDecayTrackIndexes_;
|
| 53 |
|
| 54 |
void mapTrackIdToDecayTrackIndexes();
|
| 55 |
|
| 56 |
void createTunedCollections();
|
| 57 |
|
| 58 |
void addSimulatedHistory(SimTrack const &, int parentSimTrackIndex = -1);
|
| 59 |
|
| 60 |
bool selectSimulatedTrack(SimTrack const &);
|
| 61 |
|
| 62 |
HepMC::GenParticle * getGeneratorParticle(SimTrack const &);
|
| 63 |
};
|
| 64 |
|
| 65 |
|
| 66 |
std::string HiEventTuner::MessageCategory_ = "HiEventTuner";
|
| 67 |
|
| 68 |
|
| 69 |
HiEventTuner::HiEventTuner(const edm::ParameterSet& config)
|
| 70 |
{
|
| 71 |
// Initialize global parameters
|
| 72 |
useMultipleHepMCLabels_ = config.getParameter<bool>("useMultipleHepMCLabels");
|
| 73 |
hepMCLabels_ = config.getParameter<std::vector<std::string> >("hepMCLabels");
|
| 74 |
simHitLabel_ = config.getParameter<std::string>("simHitLabel");
|
| 75 |
|
| 76 |
// Declaration of the products
|
| 77 |
produces<edm::SimTrackContainer>();
|
| 78 |
produces<edm::SimVertexContainer>();
|
| 79 |
}
|
| 80 |
|
| 81 |
|
| 82 |
void HiEventTuner::produce(edm::Event& event, const edm::EventSetup& setup)
|
| 83 |
{
|
| 84 |
// Clean the list of hepmc products
|
| 85 |
hepMCProducts_.clear();
|
| 86 |
|
| 87 |
// Collect all the HepMCProducts
|
| 88 |
edm::Handle<edm::HepMCProduct> hepMCHandle;
|
| 89 |
|
| 90 |
for (std::vector<std::string>::const_iterator source = hepMCLabels_.begin(); source != hepMCLabels_.end(); ++source)
|
| 91 |
{
|
| 92 |
if ( event.getByLabel(*source, hepMCHandle) )
|
| 93 |
{
|
| 94 |
hepMCProducts_.push_back(hepMCHandle);
|
| 95 |
edm::LogInfo (MessageCategory_) << "Using HepMC source " << *source;
|
| 96 |
if (!useMultipleHepMCLabels_) break;
|
| 97 |
}
|
| 98 |
}
|
| 99 |
|
| 100 |
// Warning and checking preconditions
|
| 101 |
if ( hepMCProducts_.empty() )
|
| 102 |
{
|
| 103 |
edm::LogWarning (MessageCategory_) << "No HepMC source found";
|
| 104 |
return;
|
| 105 |
}
|
| 106 |
else if ( hepMCProducts_.size() > 1 || useMultipleHepMCLabels_ )
|
| 107 |
{
|
| 108 |
edm::LogInfo (MessageCategory_) << "You are using more than one HepMC source.";
|
| 109 |
edm::LogInfo (MessageCategory_) << "If the labels are not in the same order as the events in the crossing frame (i.e. signal, pileup(s) ) ";
|
| 110 |
edm::LogInfo (MessageCategory_) << "or there are fewer labels than events in the crossing frame";
|
| 111 |
edm::LogInfo (MessageCategory_) << MessageCategory_ << " may try to access data in the wrong HepMCProduct and crash.";
|
| 112 |
}
|
| 113 |
|
| 114 |
// Collect all the SimTracks
|
| 115 |
event.getByLabel("mix", simHitLabel_, simTracks_);
|
| 116 |
|
| 117 |
// Collect all the simvertex from the crossing frame
|
| 118 |
event.getByLabel("mix", simHitLabel_, simVertexes_);
|
| 119 |
|
| 120 |
// Create collections of things we will put in event
|
| 121 |
tunedSimTracks_ = std::auto_ptr<edm::SimTrackContainer>( new edm::SimTrackContainer );
|
| 122 |
tunedSimVertexes_ = std::auto_ptr<edm::SimVertexContainer>( new edm::SimVertexContainer );
|
| 123 |
|
| 124 |
// Create a map between trackId and vertex (decay) index
|
| 125 |
mapTrackIdToDecayTrackIndexes();
|
| 126 |
|
| 127 |
// Create the actual tuned collection
|
| 128 |
createTunedCollections();
|
| 129 |
|
| 130 |
// Put TrackingParticles and TrackingVertices in event
|
| 131 |
event.put(tunedSimTracks_);
|
| 132 |
event.put(tunedSimVertexes_);
|
| 133 |
}
|
| 134 |
|
| 135 |
|
| 136 |
void HiEventTuner::mapTrackIdToDecayTrackIndexes()
|
| 137 |
{
|
| 138 |
std::size_t index = 0;
|
| 139 |
|
| 140 |
// Clear the association map
|
| 141 |
trackIndexToDecayTrackIndexes_.clear();
|
| 142 |
|
| 143 |
// Loop over the track collection
|
| 144 |
for (edm::SimTrackContainer::const_iterator simTrack = simTracks_->begin(); simTrack != simTracks_->end(); ++simTrack, ++index)
|
| 145 |
{
|
| 146 |
// Check for a source vertex
|
| 147 |
if (simTrack->noVertex()) continue;
|
| 148 |
|
| 149 |
// Get the source vertex
|
| 150 |
SimVertex const & sourceVertex = simVertexes_->at( simTrack->vertIndex() );
|
| 151 |
|
| 152 |
// Check if the source vertex has a parent track
|
| 153 |
if (sourceVertex.noParent()) continue;
|
| 154 |
|
| 155 |
// Get the source vertex
|
| 156 |
SimTrack const & parentTrack = simTracks_->at( sourceVertex.parentIndex() );
|
| 157 |
|
| 158 |
// Making the association between parentSimTrack and decay index
|
| 159 |
trackIndexToDecayTrackIndexes_.insert( std::make_pair(parentTrack.trackId(), index) );
|
| 160 |
}
|
| 161 |
}
|
| 162 |
|
| 163 |
|
| 164 |
void HiEventTuner::createTunedCollections()
|
| 165 |
{
|
| 166 |
// Loop over the tracks looking for the seeding tracks
|
| 167 |
for (edm::SimTrackContainer::const_iterator simTrack = simTracks_->begin(); simTrack != simTracks_->end(); ++simTrack)
|
| 168 |
if ( selectSimulatedTrack(*simTrack) ) addSimulatedHistory(*simTrack);
|
| 169 |
}
|
| 170 |
|
| 171 |
|
| 172 |
void HiEventTuner::addSimulatedHistory(SimTrack const & simTrack, int parentSimTrackIndex)
|
| 173 |
{
|
| 174 |
// Add the simTrack to the tuned collection
|
| 175 |
tunedSimTracks_->push_back(simTrack);
|
| 176 |
|
| 177 |
// Reference to the tunedSimTrack
|
| 178 |
SimTrack & tunedSimTrack = tunedSimTracks_->back();
|
| 179 |
|
| 180 |
// Add the origin vertex of a selected track
|
| 181 |
if ( !simTrack.noVertex() )
|
| 182 |
{
|
| 183 |
// Add the parent simVertex to the tunedSimVertexes collection
|
| 184 |
tunedSimVertexes_->push_back(
|
| 185 |
SimVertex(simVertexes_->at(simTrack.vertIndex()), parentSimTrackIndex)
|
| 186 |
);
|
| 187 |
|
| 188 |
// Update the tunedSimTrack vertex index
|
| 189 |
tunedSimTrack.setVertexIndex((int)tunedSimVertexes_->size());
|
| 190 |
}
|
| 191 |
|
| 192 |
// Define trackId
|
| 193 |
std::size_t simTrackIndex = simTrack.trackId();
|
| 194 |
|
| 195 |
// Loop over the decay tracks
|
| 196 |
for (
|
| 197 |
IndexToIndexes::const_iterator decaySimTrackIndex = trackIndexToDecayTrackIndexes_.lower_bound(simTrackIndex);
|
| 198 |
decaySimTrackIndex != trackIndexToDecayTrackIndexes_.upper_bound(simTrackIndex);
|
| 199 |
++decaySimTrackIndex
|
| 200 |
)
|
| 201 |
addSimulatedHistory(simTracks_->at(decaySimTrackIndex->second), (int)tunedSimTracks_->size());
|
| 202 |
}
|
| 203 |
|
| 204 |
|
| 205 |
bool HiEventTuner::selectSimulatedTrack(SimTrack const & simTrack)
|
| 206 |
{
|
| 207 |
// This is the basic rule, select SimTracks associate with a GenParticle
|
| 208 |
if (getGeneratorParticle(simTrack)) return true;
|
| 209 |
return false;
|
| 210 |
}
|
| 211 |
|
| 212 |
|
| 213 |
HepMC::GenParticle * HiEventTuner::getGeneratorParticle(SimTrack const & simTrack)
|
| 214 |
{
|
| 215 |
// HepMC holder
|
| 216 |
edm::Handle<edm::HepMCProduct> hepmc;
|
| 217 |
|
| 218 |
// Index to the generator particle
|
| 219 |
int genParticleIndex = simTrack.genpartIndex();
|
| 220 |
|
| 221 |
// Check in case there is a GenParticle associated to SimTrack
|
| 222 |
if (genParticleIndex >= 0)
|
| 223 |
{
|
| 224 |
// Get the correct HepMC product (for signal and pileup if any)
|
| 225 |
hepmc = (useMultipleHepMCLabels_) ? hepMCProducts_.at(simTrack.eventId().rawId()) : hepmc = hepMCProducts_.at(0);
|
| 226 |
|
| 227 |
// Return the pointer to the GenParticle if there is any
|
| 228 |
return hepmc->GetEvent()->barcode_to_particle(genParticleIndex);
|
| 229 |
}
|
| 230 |
|
| 231 |
return 0;
|
| 232 |
}
|
| 233 |
|
| 234 |
|
| 235 |
//define this as a plug-in
|
| 236 |
DEFINE_FWK_MODULE(HiEventTuner);
|
