1 |
bazterra |
1.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);
|