CMS 3D CMS Logo

HTTTopJetProducer.cc
Go to the documentation of this file.
4 #include "HTTTopJetProducer.h"
5 
6 using namespace edm;
7 using namespace cms;
8 using namespace reco;
9 using namespace std;
10 
11 HTTTopJetProducer::HTTTopJetProducer(edm::ParameterSet const& conf):
12  FastjetJetProducer( conf ),
13  optimalR_(false),
14  qJets_(false),
15  minFatjetPt_(200.),
16  minSubjetPt_(20.),
17  minCandPt_(200.),
18  maxFatjetAbsEta_(2.5),
19  subjetMass_(30.),
20  muCut_(0.8),
21  filtR_(0.3),
22  filtN_(5),
23  mode_(0),
24  minCandMass_(150.),
25  maxCandMass_(200.),
26  massRatioWidth_(15),
27  minM23Cut_(0.35),
28  minM13Cut_(0.2),
29  maxM13Cut_(1.3),
30  maxR_(1.5),
31  minR_(0.5),
32  rejectMinR_(false),
33  verbose_(false )
34 {
35 
36  // Read in all the options from the configuration
37  optimalR_ = conf.getParameter<bool>("optimalR");
38  qJets_ = conf.getParameter<bool>("qJets");
39  minFatjetPt_ = conf.getParameter<double>("minFatjetPt");
40  minSubjetPt_ = conf.getParameter<double>("minSubjetPt");
41  minCandPt_ = conf.getParameter<double>("minCandPt");
42  maxFatjetAbsEta_ = conf.getParameter<double>("maxFatjetAbsEta");
43  subjetMass_ = conf.getParameter<double>("subjetMass");
44  muCut_ = conf.getParameter<double>("muCut");
45  filtR_ = conf.getParameter<double>("filtR");
46  filtN_ = conf.getParameter<int>("filtN");
47  mode_ = conf.getParameter<int>("mode");
48  minCandMass_ = conf.getParameter<double>("minCandMass");
49  maxCandMass_ = conf.getParameter<double>("maxCandMass");
50  massRatioWidth_ = conf.getParameter<double>("massRatioWidth");
51  minM23Cut_ = conf.getParameter<double>("minM23Cut");
52  minM13Cut_ = conf.getParameter<double>("minM13Cut");
53  maxM13Cut_ = conf.getParameter<double>("maxM13Cut");
54  maxR_ = conf.getParameter<double>("maxR");
55  minR_ = conf.getParameter<double>("minR");
56  rejectMinR_ = conf.getParameter<bool>("rejectMinR");
57  verbose_ = conf.getParameter<bool>("verbose");
58 
59  // Create the tagger-wrapper
60  produces<HTTTopJetTagInfoCollection>();
61 
62  // Signal to the VirtualJetProducer that we have to add HTT information
64 
65  fjHEPTopTagger_ = std::auto_ptr<fastjet::HEPTopTaggerV2>(new fastjet::HEPTopTaggerV2(
66  optimalR_,
67  qJets_,
68  minSubjetPt_,
69  minCandPt_,
70  subjetMass_,
71  muCut_,
72  filtR_,
73  filtN_,
74  mode_,
75  minCandMass_,
76  maxCandMass_,
77  massRatioWidth_,
78  minM23Cut_,
79  minM13Cut_,
80  maxM13Cut_,
81  rejectMinR_
82  ));
83 
84 }
85 
86 
87 
88 
89 
91 {
92 
93  if (qJets_){
95  CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
96  fjHEPTopTagger_->set_rng(engine);
97  }
98 
100 }
101 
103 {
104 
105  if ( !doAreaFastjet_ && !doRhoFastjet_) {
106  fjClusterSeq_ = ClusterSequencePtr( new fastjet::ClusterSequence( fjInputs_, *fjJetDefinition_ ) );
107  } else if (voronoiRfact_ <= 0) {
108  fjClusterSeq_ = ClusterSequencePtr( new fastjet::ClusterSequenceArea( fjInputs_, *fjJetDefinition_ , *fjAreaDefinition_ ) );
109  } else {
110  fjClusterSeq_ = ClusterSequencePtr( new fastjet::ClusterSequenceVoronoiArea( fjInputs_, *fjJetDefinition_ , fastjet::VoronoiAreaSpec(voronoiRfact_) ) );
111  }
112 
113  //Run the jet clustering
114  vector<fastjet::PseudoJet> inclusiveJets = fjClusterSeq_->inclusive_jets(minFatjetPt_);
115 
116  if ( verbose_ ) cout << "Getting central jets" << endl;
117  // Find the transient central jets
118  vector<fastjet::PseudoJet> centralJets;
119  for (unsigned int i = 0; i < inclusiveJets.size(); i++) {
120 
121  if (inclusiveJets[i].perp() > minFatjetPt_ && fabs(inclusiveJets[i].rapidity()) < maxFatjetAbsEta_) {
122  centralJets.push_back(inclusiveJets[i]);
123  }
124  }
125 
126  fastjet::HEPTopTaggerV2 & HEPTagger = *fjHEPTopTagger_;
127 
128  vector<fastjet::PseudoJet>::iterator jetIt = centralJets.begin(), centralJetsEnd = centralJets.end();
129  if ( verbose_ )cout<<"Loop over jets"<<endl;
130  for ( ; jetIt != centralJetsEnd; ++jetIt ) {
131 
132  if (verbose_) cout << "CMS FJ jet pt: " << (*jetIt).perp() << endl;
133 
134  fastjet::PseudoJet taggedJet;
135 
136  taggedJet = HEPTagger.result(*jetIt);
137 
138  if (taggedJet != 0){
139  fjJets_.push_back(taggedJet);
140  }
141  }
142 
143 }
144 
146  const edm::EventSetup& iSetup,
148 
149 
150  // Set up output list
151  auto tagInfos = std::make_unique<HTTTopJetTagInfoCollection>();
152 
153  // Loop over jets
154  for (size_t ij=0; ij != fjJets_.size(); ij++){
155 
156  HTTTopJetProperties properties;
157  HTTTopJetTagInfo tagInfo;
158 
159  // Black magic:
160  // In the standard CA treatment the RefToBase is made from the handle directly
161  // Since we only have a OrphanHandle (the JetCollection is created by this process)
162  // we have to take the detour via the Ref
164  edm::RefToBase<reco::Jet> rtb(ref);
165 
166  fastjet::HEPTopTaggerV2Structure *s = (fastjet::HEPTopTaggerV2Structure*) fjJets_[ij].structure_non_const_ptr();
167 
168  properties.fjMass = s->fj_mass();
169  properties.fjPt = s->fj_pt();
170  properties.fjEta = s->fj_eta();
171  properties.fjPhi = s->fj_phi();
172 
173  properties.topMass = s->top_mass();
174  properties.unfilteredMass = s->unfiltered_mass();
175  properties.prunedMass = s->pruned_mass();
176  properties.fRec = s->fRec();
177  properties.massRatioPassed = s->mass_ratio_passed();
178 
179  properties.ropt = s->ropt();
180  properties.roptCalc = s->roptCalc();
181  properties.ptForRoptCalc = s->ptForRoptCalc();
182 
183  properties.tau1Unfiltered = s->Tau1Unfiltered();
184  properties.tau2Unfiltered = s->Tau2Unfiltered();
185  properties.tau3Unfiltered = s->Tau3Unfiltered();
186  properties.tau1Filtered = s->Tau1Filtered();
187  properties.tau2Filtered = s->Tau2Filtered();
188  properties.tau3Filtered = s->Tau3Filtered();
189  properties.qWeight = s->qweight();
190  properties.qEpsilon = s->qepsilon();
191  properties.qSigmaM = s->qsigmaM();
192 
193  tagInfo.insert(rtb, properties );
194  tagInfos->push_back( tagInfo );
195  }
196 
197  iEvent.put(std::move(tagInfos));
198 
199 };
200 
202 
204 
205  desc.add<bool> ("optimalR", true);
206  desc.add<bool> ("qJets", false);
207  desc.add<double>("minFatjetPt", 200.);
208  desc.add<double>("minSubjetPt", 0.);
209  desc.add<double>("minCandPt", 0.);
210  desc.add<double>("maxFatjetAbsEta", 99.);
211  desc.add<double>("subjetMass", 30.);
212  desc.add<double>("filtR", 0.3);
213  desc.add<int> ("filtN", 5);
214  desc.add<int> ("mode", 4);
215  desc.add<double>("minCandMass", 0.);
216  desc.add<double>("maxCandMass", 999999.);
217  desc.add<double>("massRatioWidth",999999.);
218  desc.add<double>("minM23Cut", 0.);
219  desc.add<double>("minM13Cut", 0.);
220  desc.add<double>("maxM13Cut", 999999.);
221  desc.add<double>("maxR", 1.5);
222  desc.add<double>("minR", 0.5);
223  desc.add<bool> ("rejectMinR", false);
224  desc.add<bool> ("verbose", false);
225 
226  desc.add<int> ("algorithm", 1); // where is it needed?
227 
231  desc.add<std::string> ("jetCollInstanceName","SubJets");
233 
234  descriptions.add("HTTTopJetProducer",desc);
235 
236 }
237 
238 //define this as a plug-in
T getParameter(std::string const &) const
std::auto_ptr< fastjet::HEPTopTaggerV2 > fjHEPTopTagger_
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:127
std::vector< fastjet::PseudoJet > fjJets_
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
virtual CLHEP::HepRandomEngine & getEngine(StreamID const &)=0
Use this engine in event methods.
void insert(const edm::RefToBase< Jet > &jet, const HTTTopJetProperties &properties)
void produce(edm::Event &iEvent, const edm::EventSetup &iSetup) override
std::vector< fastjet::PseudoJet > fjInputs_
int iEvent
Definition: GenABIO.cc:230
ClusterSequencePtr fjClusterSeq_
void produce(edm::Event &iEvent, const edm::EventSetup &iSetup) override
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
ParameterDescriptionBase * add(U const &iLabel, T const &value)
void addHTTTopJetTagInfoCollection(edm::Event &iEvent, const edm::EventSetup &iSetup, edm::OrphanHandle< reco::BasicJetCollection > &oh) override
static void fillDescriptionsFromFastJetProducer(edm::ParameterSetDescription &desc)
void add(std::string const &label, ParameterSetDescription const &psetDescription)
T perp() const
Magnitude of transverse component.
fixed size matrix
HLT enums.
StreamID streamID() const
Definition: Event.h:86
AreaDefinitionPtr fjAreaDefinition_
static void fillDescriptionsFromVirtualJetProducer(edm::ParameterSetDescription &desc)
def move(src, dest)
Definition: eostools.py:510
void runAlgorithm(edm::Event &iEvent, const edm::EventSetup &iSetup) override
boost::shared_ptr< fastjet::ClusterSequence > ClusterSequencePtr