CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
FastjetJetProducer.cc
Go to the documentation of this file.
1 //
3 // FastjetJetProducer
4 // ------------------
5 //
6 // 04/21/2009 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
8 
10 
19 
25 
28 
29 #include "fastjet/SISConePlugin.hh"
30 #include "fastjet/CMSIterativeConePlugin.hh"
31 #include "fastjet/ATLASConePlugin.hh"
32 #include "fastjet/CDFMidPointPlugin.hh"
33 #include "fastjet/tools/Filter.hh"
34 #include "fastjet/tools/Pruner.hh"
35 #include "fastjet/tools/MassDropTagger.hh"
37 
38 #include <iostream>
39 #include <memory>
40 #include <algorithm>
41 #include <limits>
42 #include <cmath>
43 //#include <fstream>
44 
45 using namespace std;
46 
47 
48 
50 // construction / destruction
52 
53 //______________________________________________________________________________
55  : VirtualJetProducer( iConfig ),
56  useMassDropTagger_(false),
57  useFiltering_(false),
58  useTrimming_(false),
59  usePruning_(false),
60  useCMSBoostedTauSeedingAlgorithm_(false),
61  muCut_(-1.0),
62  yCut_(-1.0),
63  rFilt_(-1.0),
64  nFilt_(-1),
65  trimPtFracMin_(-1.0),
66  zCut_(-1.0),
67  RcutFactor_(-1.0)
68 {
69 
70  if ( iConfig.exists("UseOnlyVertexTracks") )
71  useOnlyVertexTracks_ = iConfig.getParameter<bool>("UseOnlyVertexTracks");
72  else
73  useOnlyVertexTracks_ = false;
74 
75  if ( iConfig.exists("UseOnlyOnePV") )
76  useOnlyOnePV_ = iConfig.getParameter<bool>("UseOnlyOnePV");
77  else
78  useOnlyOnePV_ = false;
79 
80  if ( iConfig.exists("DzTrVtxMax") )
81  dzTrVtxMax_ = iConfig.getParameter<double>("DzTrVtxMax");
82  else
83  dzTrVtxMax_ = 999999.;
84  if ( iConfig.exists("DxyTrVtxMax") )
85  dxyTrVtxMax_ = iConfig.getParameter<double>("DxyTrVtxMax");
86  else
87  dxyTrVtxMax_ = 999999.;
88  if ( iConfig.exists("MinVtxNdof") )
89  minVtxNdof_ = iConfig.getParameter<int>("MinVtxNdof");
90  else
91  minVtxNdof_ = 5;
92  if ( iConfig.exists("MaxVtxZ") )
93  maxVtxZ_ = iConfig.getParameter<double>("MaxVtxZ");
94  else
95  maxVtxZ_ = 15;
96 
97 
98  if ( iConfig.exists("useFiltering") ||
99  iConfig.exists("useTrimming") ||
100  iConfig.exists("usePruning") ||
101  iConfig.exists("useMassDropTagger") ||
102  iConfig.exists("useCMSBoostedTauSeedingAlgorithm") ) {
103  useMassDropTagger_=false;
104  useFiltering_=false;
105  useTrimming_=false;
106  usePruning_=false;
108  rFilt_=-1.0;
109  nFilt_=-1;
110  trimPtFracMin_=-1.0;
111  zCut_=-1.0;
112  RcutFactor_=-1.0;
113  muCut_=-1.0;
114  yCut_=-1.0;
115  subjetPtMin_ = -1.0;
116  muMin_ = -1.0;
117  muMax_ = -1.0;
118  yMin_ = -1.0;
119  yMax_ = -1.0;
120  dRMin_ = -1.0;
121  dRMax_ = -1.0;
122  maxDepth_ = -1;
123  useExplicitGhosts_ = true;
124 
125  if ( iConfig.exists("useMassDropTagger") ) {
126  useMassDropTagger_ = iConfig.getParameter<bool>("useMassDropTagger");
127  muCut_ = iConfig.getParameter<double>("muCut");
128  yCut_ = iConfig.getParameter<double>("yCut");
129  }
130 
131  if ( iConfig.exists("useFiltering") ) {
132  useFiltering_ = iConfig.getParameter<bool>("useFiltering");
133  rFilt_ = iConfig.getParameter<double>("rFilt");
134  nFilt_ = iConfig.getParameter<int>("nFilt");
135  }
136 
137  if ( iConfig.exists("useTrimming") ) {
138  useTrimming_ = iConfig.getParameter<bool>("useTrimming");
139  rFilt_ = iConfig.getParameter<double>("rFilt");
140  trimPtFracMin_ = iConfig.getParameter<double>("trimPtFracMin");
141  }
142 
143  if ( iConfig.exists("usePruning") ) {
144  usePruning_ = iConfig.getParameter<bool>("usePruning");
145  zCut_ = iConfig.getParameter<double>("zcut");
146  RcutFactor_ = iConfig.getParameter<double>("rcut_factor");
147  nFilt_ = iConfig.getParameter<int>("nFilt");
148  }
149 
150  if ( iConfig.exists("useCMSBoostedTauSeedingAlgorithm") ) {
151  useCMSBoostedTauSeedingAlgorithm_ = iConfig.getParameter<bool>("useCMSBoostedTauSeedingAlgorithm");
152  subjetPtMin_ = iConfig.getParameter<double>("subjetPtMin");
153  muMin_ = iConfig.getParameter<double>("muMin");
154  muMax_ = iConfig.getParameter<double>("muMax");
155  yMin_ = iConfig.getParameter<double>("yMin");
156  yMax_ = iConfig.getParameter<double>("yMax");
157  dRMin_ = iConfig.getParameter<double>("dRMin");
158  dRMax_ = iConfig.getParameter<double>("dRMax");
159  maxDepth_ = iConfig.getParameter<int>("maxDepth");
160  }
161 
162  }
163 
164  input_chrefcand_token_ = consumes<edm::View<reco::RecoChargedRefCandidate> >(src_);
165 
166 }
167 
168 
169 //______________________________________________________________________________
171 {
172 }
173 
174 
176 // implementation of member functions
178 
180 {
181 
182  // for everything but track jets
183  if (!makeTrackJet(jetTypeE)) {
184 
185  // use the default production from one collection
186  VirtualJetProducer::produce( iEvent, iSetup );
187 
188  } else { // produce trackjets from tracks grouped per primary vertex
189 
190  produceTrackJets(iEvent, iSetup);
191 
192  }
193 
194 }
195 
196 
198 {
199 
200  // read in the track candidates
202  iEvent.getByToken(input_chrefcand_token_, inputsHandle);
203 
204  // make collection with pointers so we can play around with it
205  std::vector<edm::Ptr<reco::RecoChargedRefCandidate> > allInputs;
206  std::vector<edm::Ptr<reco::Candidate> > origInputs;
207  for (size_t i = 0; i < inputsHandle->size(); ++i) {
208  allInputs.push_back(inputsHandle->ptrAt(i));
209  origInputs.push_back(inputsHandle->ptrAt(i));
210  }
211 
212  // read in the PV collection
214  iEvent.getByToken(input_vertex_token_, pvCollection);
215  // define the overall output jet container
216  std::auto_ptr<std::vector<reco::TrackJet> > jets(new std::vector<reco::TrackJet>() );
217 
218  // loop over the good vertices, clustering for each vertex separately
219  for (reco::VertexCollection::const_iterator itVtx = pvCollection->begin(); itVtx != pvCollection->end(); ++itVtx) {
220  if (itVtx->isFake() || itVtx->ndof() < minVtxNdof_ || fabs(itVtx->z()) > maxVtxZ_) continue;
221 
222  // clear the intermediate containers
223  inputs_.clear();
224  fjInputs_.clear();
225  fjJets_.clear();
226 
227  // if only vertex-associated tracks should be used
228  if (useOnlyVertexTracks_) {
229  // loop over the tracks associated to the vertex
230  for (reco::Vertex::trackRef_iterator itTr = itVtx->tracks_begin(); itTr != itVtx->tracks_end(); ++itTr) {
231  // whether a match was found in the track candidate input
232  bool found = false;
233  // loop over input track candidates
234  for (std::vector<edm::Ptr<reco::RecoChargedRefCandidate> >::iterator itIn = allInputs.begin(); itIn != allInputs.end(); ++itIn) {
235  // match the input track candidate to the track from the vertex
236  reco::TrackRef trref(itTr->castTo<reco::TrackRef>());
237  // check if the tracks match
238  if ((*itIn)->track() == trref) {
239  found = true;
240  // add this track candidate to the input for clustering
241  inputs_.push_back(*itIn);
242  // erase the track candidate from the total list of input, so we don't reuse it later
243  allInputs.erase(itIn);
244  // found the candidate track corresponding to the vertex track, so stop the loop
245  break;
246  } // end if match found
247  } // end loop over input tracks
248  // give an info message in case no match is found (can happen if candidates are subset of tracks used for clustering)
249  if (!found) edm::LogInfo("FastjetTrackJetProducer") << "Ignoring a track at vertex which is not in input track collection!";
250  } // end loop over tracks associated to vertex
251  // if all inpt track candidates should be used
252  } else {
253  // loop over input track candidates
254  for (std::vector<edm::Ptr<reco::RecoChargedRefCandidate> >::iterator itIn = allInputs.begin(); itIn != allInputs.end(); ++itIn) {
255  // check if the track is close enough to the vertex
256  float dz = (*itIn)->track()->dz(itVtx->position());
257  float dxy = (*itIn)->track()->dxy(itVtx->position());
258  if (fabs(dz) > dzTrVtxMax_) continue;
259  if (fabs(dxy) > dxyTrVtxMax_) continue;
260  bool closervtx = false;
261  // now loop over the good vertices a second time
262  for (reco::VertexCollection::const_iterator itVtx2 = pvCollection->begin(); itVtx2 != pvCollection->end(); ++itVtx2) {
263  if (itVtx->isFake() || itVtx->ndof() < minVtxNdof_ || fabs(itVtx->z()) > maxVtxZ_) continue;
264  // and check this track is closer to any other vertex (if more than 1 vertex considered)
265  if (!useOnlyOnePV_ &&
266  itVtx != itVtx2 &&
267  fabs((*itIn)->track()->dz(itVtx2->position())) < fabs(dz)) {
268  closervtx = true;
269  break; // 1 closer vertex makes the track already not matched, so break
270  }
271  }
272  // don't add this track if another vertex is found closer
273  if (closervtx) continue;
274  // add this track candidate to the input for clustering
275  inputs_.push_back(*itIn);
276  // erase the track candidate from the total list of input, so we don't reuse it later
277  allInputs.erase(itIn);
278  // take a step back in the loop since we just erased
279  --itIn;
280  }
281  }
282 
283  // convert candidates in inputs_ to fastjet::PseudoJets in fjInputs_
284  fjInputs_.reserve(inputs_.size());
285  inputTowers();
286  LogDebug("FastjetTrackJetProducer") << "Inputted towers\n";
287 
288  // run algorithm, using fjInputs_, modifying fjJets_ and allocating fjClusterSeq_
289  runAlgorithm(iEvent, iSetup);
290  LogDebug("FastjetTrackJetProducer") << "Ran algorithm\n";
291 
292  // convert our jets and add to the overall jet vector
293  for (unsigned int ijet=0;ijet<fjJets_.size();++ijet) {
294  // get the constituents from fastjet
295  std::vector<fastjet::PseudoJet> fjConstituents = sorted_by_pt(fjClusterSeq_->constituents(fjJets_[ijet]));
296  // convert them to CandidatePtr vector
297  std::vector<reco::CandidatePtr> constituents = getConstituents(fjConstituents);
298  // fill the trackjet
300  // write the specifics to the jet (simultaneously sets 4-vector, vertex).
301  writeSpecific( jet,
302  reco::Particle::LorentzVector(fjJets_[ijet].px(), fjJets_[ijet].py(), fjJets_[ijet].pz(), fjJets_[ijet].E()),
303  vertex_, constituents, iSetup);
304  jet.setJetArea(0);
305  jet.setPileup(0);
306  jet.setPrimaryVertex(edm::Ref<reco::VertexCollection>(pvCollection, (int) (itVtx-pvCollection->begin())));
307  jet.setVertex(itVtx->position());
308  jets->push_back(jet);
309  }
310 
311  if (useOnlyOnePV_) break; // stop vertex loop if only one vertex asked for
312  } // end loop over vertices
313 
314  // put the jets in the collection
315  LogDebug("FastjetTrackJetProducer") << "Put " << jets->size() << " jets in the event.\n";
316  iEvent.put(jets);
317 
318 }
319 
320 
321 //______________________________________________________________________________
323 {
324  // run algorithm
325  /*
326  fjInputs_.clear();
327  double px, py , pz, E;
328  string line;
329  std::ifstream fin("dump3.txt");
330  while (getline(fin, line)){
331  if (line == "#END") break;
332  if (line.substr(0,1) == "#") {continue;}
333  istringstream istr(line);
334  istr >> px >> py >> pz >> E;
335  // create a fastjet::PseudoJet with these components and put it onto
336  // back of the input_particles vector
337  fastjet::PseudoJet j(px,py,pz,E);
338  //if ( fabs(j.rap()) < inputEtaMax )
339  fjInputs_.push_back(fastjet::PseudoJet(px,py,pz,E));
340  }
341  fin.close();
342  */
343 
344  if ( !doAreaFastjet_ && !doRhoFastjet_) {
345  fjClusterSeq_ = ClusterSequencePtr( new fastjet::ClusterSequence( fjInputs_, *fjJetDefinition_ ) );
346  } else if (voronoiRfact_ <= 0) {
347  fjClusterSeq_ = ClusterSequencePtr( new fastjet::ClusterSequenceArea( fjInputs_, *fjJetDefinition_ , *fjAreaDefinition_ ) );
348  } else {
349  fjClusterSeq_ = ClusterSequencePtr( new fastjet::ClusterSequenceVoronoiArea( fjInputs_, *fjJetDefinition_ , fastjet::VoronoiAreaSpec(voronoiRfact_) ) );
350  }
351 
353  fjJets_ = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(jetPtMin_));
354  } else {
355  fjJets_.clear();
356  std::vector<fastjet::PseudoJet> tempJets = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(jetPtMin_));
357 
358  fastjet::MassDropTagger md_tagger( muCut_, yCut_ );
359  fastjet::contrib::CMSBoostedTauSeedingAlgorithm tau_tagger( subjetPtMin_, muMin_, muMax_, yMin_, yMax_, dRMin_, dRMax_, maxDepth_, verbosity_ );
360  fastjet::Filter trimmer( fastjet::Filter(fastjet::JetDefinition(fastjet::kt_algorithm, rFilt_), fastjet::SelectorPtFractionMin(trimPtFracMin_)));
361  fastjet::Filter filter( fastjet::Filter(fastjet::JetDefinition(fastjet::cambridge_algorithm, rFilt_), fastjet::SelectorNHardest(nFilt_)));
362  fastjet::Pruner pruner(fastjet::cambridge_algorithm, zCut_, RcutFactor_);
363 
364  std::vector<fastjet::Transformer const *> transformers;
365 
366  if ( useMassDropTagger_ ) {
367  transformers.push_back(&md_tagger);
368  }
370  transformers.push_back(&tau_tagger);
371  }
372  if ( useTrimming_ ) {
373  transformers.push_back(&trimmer);
374  }
375  if ( useFiltering_ ) {
376  transformers.push_back(&filter);
377  }
378  if ( usePruning_ ) {
379  transformers.push_back(&pruner);
380  }
381 
382  for ( std::vector<fastjet::PseudoJet>::const_iterator ijet = tempJets.begin(),
383  ijetEnd = tempJets.end(); ijet != ijetEnd; ++ijet ) {
384 
385  fastjet::PseudoJet transformedJet = *ijet;
386  bool passed = true;
387  for ( std::vector<fastjet::Transformer const *>::const_iterator itransf = transformers.begin(),
388  itransfEnd = transformers.end(); itransf != itransfEnd; ++itransf ) {
389  if ( transformedJet != 0 ) {
390  transformedJet = (**itransf)(transformedJet);
391  } else {
392  passed=false;
393  }
394  }
395  if ( passed ) {
396  fjJets_.push_back( transformedJet );
397  }
398  }
399  }
400 
401 }
402 
403 
404 
406 // define as cmssw plugin
408 
410 
#define LogDebug(id)
double muMin_
for CMSBoostedTauSeedingAlgorithm : subjet pt min
JetType::Type jetTypeE
T getParameter(std::string const &) const
double yMin_
for CMSBoostedTauSeedingAlgorithm : max mass-drop
bool useFiltering_
Mass-drop tagging for boosted Higgs.
int i
Definition: DBlmapReader.cc:9
double rFilt_
for mass-drop tagging, symmetry cut: min(pt1^2,pt2^2) * dR(1,2) / mjet &gt; ycut
tuple Filter
Definition: Filter_cff.py:5
reco::Particle::Point vertex_
double dRMax_
for CMSBoostedTauSeedingAlgorithm : min dR
double subjetPtMin_
for pruning: constituent dR * pt/2m &lt; rcut_factor
virtual std::vector< reco::CandidatePtr > getConstituents(const std::vector< fastjet::PseudoJet > &fjConstituents)
std::vector< fastjet::PseudoJet > fjJets_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:434
virtual void setPileup(float fEnergy)
Set pileup energy contribution as calculated by algorithm.
Definition: Jet.h:108
virtual void inputTowers()
virtual void runAlgorithm(edm::Event &iEvent, const edm::EventSetup &iSetup)
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
bool exists(std::string const &parameterName) const
checks if a parameter exists
double RcutFactor_
for pruning: constituent minimum pt fraction of parent cluster
double yMax_
for CMSBoostedTauSeedingAlgorithm : min asymmetry
FastjetJetProducer(const edm::ParameterSet &iConfig)
bool useTrimming_
Jet filtering technique.
bool makeTrackJet(const JetType::Type &fTag)
virtual void setJetArea(float fArea)
set jet area
Definition: Jet.h:103
std::vector< fastjet::PseudoJet > fjInputs_
virtual void produce(edm::Event &iEvent, const edm::EventSetup &iSetup)
double trimPtFracMin_
for filtering, pruning: number of subjets expected
int iEvent
Definition: GenABIO.cc:230
bool usePruning_
Jet trimming technique.
std::vector< edm::Ptr< reco::Candidate > > inputs_
bool useCMSBoostedTauSeedingAlgorithm_
Jet pruning technique.
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:116
ClusterSequencePtr fjClusterSeq_
virtual void produceTrackJets(edm::Event &iEvent, const edm::EventSetup &iSetup)
vector< PseudoJet > jets
edm::EDGetTokenT< edm::View< reco::RecoChargedRefCandidate > > input_chrefcand_token_
for CMSBoostedTauSeedingAlgorithm : max depth for descending into clustering sequence ...
double muMax_
for CMSBoostedTauSeedingAlgorithm : min mass-drop
double muCut_
algorithm for seeding reconstruction of boosted Taus (similar to mass-drop tagging) ...
void setPrimaryVertex(const reco::VertexRef &vtx)
set associated primary vertex
Definition: TrackJet.cc:80
double dRMin_
for CMSBoostedTauSeedingAlgorithm : max asymmetry
double zCut_
for trimming: constituent minimum pt fraction of full jet
virtual void setVertex(const Point &vertex)
set vertex
Jets made out of tracks.
Definition: TrackJet.h:27
int nFilt_
for filtering, trimming: dR scale of sub-clustering
int maxDepth_
for CMSBoostedTauSeedingAlgorithm : max dR
edm::EDGetTokenT< reco::VertexCollection > input_vertex_token_
virtual void produce(edm::Event &iEvent, const edm::EventSetup &iSetup)
double yCut_
for mass-drop tagging, m0/mjet (m0 = mass of highest mass subjet)
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector&lt;TrackRef&gt;
Definition: Vertex.h:37
AreaDefinitionPtr fjAreaDefinition_
volatile std::atomic< bool > shutdown_flag false
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Particle.h:27
boost::shared_ptr< fastjet::ClusterSequence > ClusterSequencePtr
void writeSpecific(reco::CaloJet &jet, reco::Particle::LorentzVector const &p4, reco::Particle::Point const &point, std::vector< reco::CandidatePtr > const &constituents, edm::EventSetup const &c)
Definition: JetSpecific.cc:41