CMS 3D CMS Logo
PileupJPTJetIdAlgo.cc
Go to the documentation of this file.
1 #include <fstream>
2 #include <iomanip>
3 #include <iostream>
4 
5 // DataFormat //
6 
7 #include "DataFormats/Common/interface/Handle.h"
8 #include "DataFormats/Provenance/interface/Provenance.h"
9 
10 #include "DataFormats/JetReco/interface/Jet.h"
11 #include "DataFormats/JetReco/interface/JetTracksAssociation.h"
12 
13 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
14 #include "DataFormats/JetReco/interface/CaloJet.h"
15 
16 #include "DataFormats/JetReco/interface/JPTJetCollection.h"
17 #include "DataFormats/JetReco/interface/JPTJet.h"
18 
19 #include "DataFormats/JetReco/interface/JetID.h"
20 
21 #include "DataFormats/Math/interface/deltaR.h"
22 #include "DataFormats/Math/interface/LorentzVector.h"
23 #include "DataFormats/Math/interface/Vector3D.h"
24 
25 // FWCore //
26 
27 #include "FWCore/Framework/interface/Event.h"
28 #include "FWCore/Framework/interface/MakerMacros.h"
29 #include "FWCore/Framework/interface/ESHandle.h"
30 #include "FWCore/ParameterSet/interface/ParameterSet.h"
31 
32 // TrackingTools //
33 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
34 
35 // Constants, Math //
36 #include "CLHEP/Units/GlobalPhysicalConstants.h"
37 #include "Math/GenVector/VectorUtil.h"
38 #include "Math/GenVector/PxPyPzE4D.h"
39 
40 #include "CommonTools/MVAUtils/interface/TMVAZipReader.h"
41 
42 #include <vector>
43 
45 
46 #include "RecoJets/JetProducers/interface/PileupJPTJetIdAlgo.h"
47 
49 using namespace std;
50 namespace cms {
51 
52  PileupJPTJetIdAlgo::PileupJPTJetIdAlgo(const edm::ParameterSet& iConfig) {
53  verbosity = iConfig.getParameter<int>("Verbosity");
54  tmvaWeights_ = edm::FileInPath(iConfig.getParameter<std::string>("tmvaWeightsCentral")).fullPath();
55  tmvaWeightsF_ = edm::FileInPath(iConfig.getParameter<std::string>("tmvaWeightsForward")).fullPath();
56  tmvaMethod_ = iConfig.getParameter<std::string>("tmvaMethod");
57  }
58 
59  PileupJPTJetIdAlgo::~PileupJPTJetIdAlgo() {
60  // std::cout<<" JetPlusTrack destructor "<<std::endl;
61  }
62 
63  void PileupJPTJetIdAlgo::bookMVAReader() {
64  // Read TMVA tree
65  // std::string mvaw = "RecoJets/JetProducers/data/TMVAClassification_BDTG.weights.xml";
66  // std::string mvawF = "RecoJets/JetProducers/data/TMVAClassification_BDTG.weights_F.xml";
67  // tmvaWeights_ = edm::FileInPath(mvaw).fullPath();
68  // tmvaWeightsF_ = edm::FileInPath(mvawF).fullPath();
69  // tmvaMethod_ = "BDTG method";
70 
71  if (verbosity > 0)
72  std::cout << " TMVA method " << tmvaMethod_.c_str() << " " << tmvaWeights_.c_str() << std::endl;
73  reader_ = new TMVA::Reader("!Color:!Silent");
74 
75  reader_->AddVariable("Nvtx", &Nvtx);
76  reader_->AddVariable("PtJ", &PtJ);
77  reader_->AddVariable("EtaJ", &EtaJ);
78  reader_->AddVariable("Beta", &Beta);
79  reader_->AddVariable("MultCalo", &MultCalo);
80  reader_->AddVariable("dAxis1c", &dAxis1c);
81  reader_->AddVariable("dAxis2c", &dAxis2c);
82  reader_->AddVariable("MultTr", &MultTr);
83  reader_->AddVariable("dAxis1t", &dAxis1t);
84  reader_->AddVariable("dAxis2t", &dAxis2t);
85 
86  reco::details::loadTMVAWeights(reader_, tmvaMethod_, tmvaWeights_);
87 
88  //std::cout<<" Method booked "<<std::endl;
89 
90  readerF_ = new TMVA::Reader("!Color:!Silent");
91 
92  readerF_->AddVariable("Nvtx", &Nvtx);
93  readerF_->AddVariable("PtJ", &PtJ);
94  readerF_->AddVariable("EtaJ", &EtaJ);
95  readerF_->AddVariable("MultCalo", &MultCalo);
96  readerF_->AddVariable("dAxis1c", &dAxis1c);
97  readerF_->AddVariable("dAxis2c", &dAxis2c);
98 
99  reco::details::loadTMVAWeights(readerF_, tmvaMethod_, tmvaWeightsF_);
100 
101  // std::cout<<" Method booked F "<<std::endl;
102  }
103 
104  float PileupJPTJetIdAlgo::fillJPTBlock(const reco::JPTJet* jet) {
105  if (verbosity > 0) {
106  std::cout << "================================ JET LOOP ==================== " << std::endl;
107  std::cout << "================ jetPt = " << (*jet).pt() << std::endl;
108  std::cout << "================ jetEta = " << (*jet).eta() << std::endl;
109  }
110 
111  edm::RefToBase<reco::Jet> jptjetRef = jet->getCaloJetRef();
112  reco::CaloJet const* rawcalojet = dynamic_cast<reco::CaloJet const*>(&*jptjetRef);
113 
114  int ncalotowers = 0.;
115  double sumpt = 0.;
116  double dphi2 = 0.;
117  double deta2 = 0.;
118  double dphi1 = 0.;
119  double dphideta = 0.;
120  double deta1 = 0.;
121  double EE = 0.;
122  double HE = 0.;
123  double EELong = 0.;
124  double EEShort = 0.;
125 
126  std::vector<CaloTowerPtr> calotwrs = (*rawcalojet).getCaloConstituents();
127 
128  if (verbosity > 0) {
129  std::cout << "======= CaloTowerPtr DONE == " << std::endl;
130  }
131 
132  for (std::vector<CaloTowerPtr>::const_iterator icalot = calotwrs.begin(); icalot != calotwrs.end(); icalot++) {
133  ncalotowers++;
134 
135  double deta = (*jet).eta() - (*icalot)->eta();
136  double dphi = (*jet).phi() - (*icalot)->phi();
137 
138  if (dphi > M_PI)
139  dphi = dphi - 2. * M_PI;
140  if (dphi < -1. * M_PI)
141  dphi = dphi + 2. * M_PI;
142 
143  if (verbosity > 0)
144  std::cout << " CaloTower jet eta " << (*jet).eta() << " tower eta " << (*icalot)->eta() << " jet phi "
145  << (*jet).phi() << " tower phi " << (*icalot)->phi() << " dphi " << dphi << " " << (*icalot)->pt()
146  << " ieta " << (*icalot)->ieta() << " " << abs((*icalot)->ieta()) << std::endl;
147 
148  if (abs((*icalot)->ieta()) < 30)
149  EE = EE + (*icalot)->emEnergy();
150  if (abs((*icalot)->ieta()) < 30)
151  HE = HE + (*icalot)->hadEnergy();
152  if (abs((*icalot)->ieta()) > 29)
153  EELong = EELong + (*icalot)->emEnergy();
154  if (abs((*icalot)->ieta()) > 29)
155  EEShort = EEShort + (*icalot)->hadEnergy();
156 
157  sumpt = sumpt + (*icalot)->pt();
158  dphi2 = dphi2 + dphi * dphi * (*icalot)->pt();
159  deta2 = deta2 + deta * deta * (*icalot)->pt();
160  dphi1 = dphi1 + dphi * (*icalot)->pt();
161  deta1 = deta1 + deta * (*icalot)->pt();
162  dphideta = dphideta + dphi * deta * (*icalot)->pt();
163 
164  } // calojet constituents
165 
166  if (sumpt > 0.) {
167  deta1 = deta1 / sumpt;
168  dphi1 = dphi1 / sumpt;
169  deta2 = deta2 / sumpt;
170  dphi2 = dphi2 / sumpt;
171  dphideta = dphideta / sumpt;
172  }
173 
174  // W.r.t. principal axis
175 
176  double detavar = deta2 - deta1 * deta1;
177  double dphivar = dphi2 - dphi1 * dphi1;
178  double dphidetacov = dphideta - deta1 * dphi1;
179 
180  double det = (detavar - dphivar) * (detavar - dphivar) + 4 * dphidetacov * dphidetacov;
181  det = sqrt(det);
182  double x1 = (detavar + dphivar + det) / 2.;
183  double x2 = (detavar + dphivar - det) / 2.;
184 
185  if (verbosity > 0)
186  std::cout << " ncalo " << ncalotowers << " deta2 " << deta2 << " dphi2 " << dphi2 << " deta1 " << deta1
187  << " dphi1 " << dphi1 << " detavar " << detavar << " dphivar " << dphivar << " dphidetacov "
188  << dphidetacov << " sqrt(det) " << sqrt(det) << " x1 " << x1 << " x2 " << x2 << std::endl;
189 
190  // For jets with |eta|<2 take also tracks shape
191  int ntracks = 0.;
192  double sumpttr = 0.;
193  double dphitr2 = 0.;
194  double detatr2 = 0.;
195  double dphitr1 = 0.;
196  double detatr1 = 0.;
197  double dphidetatr = 0.;
198 
199  const reco::TrackRefVector pioninin = (*jet).getPionsInVertexInCalo();
200 
201  for (reco::TrackRefVector::const_iterator it = pioninin.begin(); it != pioninin.end(); it++) {
202  if ((*it)->pt() > 0.5 && ((*it)->ptError() / (*it)->pt()) < 0.05) {
203  ntracks++;
204  sumpttr = sumpttr + (*it)->pt();
205  double deta = (*jet).eta() - (*it)->eta();
206  double dphi = (*jet).phi() - (*it)->phi();
207 
208  if (dphi > M_PI)
209  dphi = dphi - 2. * M_PI;
210  if (dphi < -1. * M_PI)
211  dphi = dphi + 2. * M_PI;
212 
213  dphitr2 = dphitr2 + dphi * dphi * (*it)->pt();
214  detatr2 = detatr2 + deta * deta * (*it)->pt();
215  dphitr1 = dphitr1 + dphi * (*it)->pt();
216  detatr1 = detatr1 + deta * (*it)->pt();
217  dphidetatr = dphidetatr + dphi * deta * (*it)->pt();
218  if (verbosity > 0)
219  std::cout << " Tracks-in-in " << (*it)->pt() << " " << (*it)->eta() << " " << (*it)->phi() << " in jet "
220  << (*jet).eta() << " " << (*jet).phi() << " jet pt " << (*jet).pt() << std::endl;
221  }
222  } // pioninin
223 
224  const reco::TrackRefVector pioninout = (*jet).getPionsInVertexOutCalo();
225 
226  for (reco::TrackRefVector::const_iterator it = pioninout.begin(); it != pioninout.end(); it++) {
227  if ((*it)->pt() > 0.5 && ((*it)->ptError() / (*it)->pt()) < 0.05) {
228  ntracks++;
229  sumpttr = sumpttr + (*it)->pt();
230  double deta = (*jet).eta() - (*it)->eta();
231  double dphi = (*jet).phi() - (*it)->phi();
232 
233  if (dphi > M_PI)
234  dphi = dphi - 2. * M_PI;
235  if (dphi < -1. * M_PI)
236  dphi = dphi + 2. * M_PI;
237 
238  dphitr2 = dphitr2 + dphi * dphi * (*it)->pt();
239  detatr2 = detatr2 + deta * deta * (*it)->pt();
240  dphitr1 = dphitr1 + dphi * (*it)->pt();
241  detatr1 = detatr1 + deta * (*it)->pt();
242  dphidetatr = dphidetatr + dphi * deta * (*it)->pt();
243  if (verbosity > 0)
244  std::cout << " Tracks-in-in " << (*it)->pt() << " " << (*it)->eta() << " " << (*it)->phi() << " in jet "
245  << (*jet).eta() << " " << (*jet).phi() << " jet pt " << (*jet).pt() << std::endl;
246  }
247  } // pioninout
248 
249  if (verbosity > 0)
250  std::cout << " Number of tracks in-in and in-out " << pioninin.size() << " " << pioninout.size() << std::endl;
251 
252  if (sumpttr > 0.) {
253  detatr1 = detatr1 / sumpttr;
254  dphitr1 = dphitr1 / sumpttr;
255  detatr2 = detatr2 / sumpttr;
256  dphitr2 = dphitr2 / sumpttr;
257  dphidetatr = dphidetatr / sumpttr;
258  }
259 
260  // W.r.t. principal axis
261 
262  double detavart = detatr2 - detatr1 * detatr1;
263  double dphivart = dphitr2 - dphitr1 * dphitr1;
264  double dphidetacovt = dphidetatr - detatr1 * dphitr1;
265 
266  double dettr = (detavart - dphivart) * (detavart - dphivart) + 4 * dphidetacovt * dphidetacovt;
267  dettr = sqrt(dettr);
268  double x1tr = (detavart + dphivart + dettr) / 2.;
269  double x2tr = (detavart + dphivart - dettr) / 2.;
270 
271  if (verbosity > 0)
272  std::cout << " ntracks " << ntracks << " detatr2 " << detatr2 << " dphitr2 " << dphitr2 << " detatr1 " << detatr1
273  << " dphitr1 " << dphitr1 << " detavart " << detavart << " dphivart " << dphivart << " dphidetacovt "
274  << dphidetacovt << " sqrt(det) " << sqrt(dettr) << " x1tr " << x1tr << " x2tr " << x2tr << std::endl;
275 
276  // Multivariate analisis
277  PtJ = (*jet).pt();
278  EtaJ = (*jet).eta();
279  Beta = (*jet).getSpecific().Zch;
280  dAxis2c = x2;
281  dAxis1c = x1;
282  dAxis2t = x2tr;
283  dAxis1t = x1tr;
284  MultCalo = ncalotowers;
285  MultTr = ntracks;
286 
287  float mva_ = 1.;
288  if (fabs(EtaJ) < 2.6) {
289  mva_ = reader_->EvaluateMVA("BDTG method");
290  // std::cout<<" MVA analysis "<<mva_<<std::endl;
291  } else {
292  mva_ = readerF_->EvaluateMVA("BDTG method");
293  // std::cout<<" MVA analysis Forward "<<mva_<<std::endl;
294  }
295  if (verbosity > 0)
296  std::cout << "======= Computed MVA = " << mva_ << std::endl;
297  return mva_;
298  } // fillJPTBlock
299 } // namespace cms
reco::JPTJet
Jets made from CaloJets corrected for ZSP and tracks.
Definition: JPTJet.h:28
HIPAlignmentAlgorithm_cfi.verbosity
verbosity
Definition: HIPAlignmentAlgorithm_cfi.py:7
JetTracksAssociation.h
CaloJetCollection.h
reco::CaloJet
Jets made from CaloTowers.
Definition: CaloJet.h:27
Handle.h
CaloJet.h
JetID.h
ESHandle.h
testProducerWithPsetDescEmpty_cfi.x2
x2
Definition: testProducerWithPsetDescEmpty_cfi.py:28
gather_cfg.cout
cout
Definition: gather_cfg.py:144
edm::RefVector::begin
const_iterator begin() const
Initialize an iterator over the RefVector.
Definition: RefVector.h:223
Jet.h
edm::RefVector< TrackCollection >
vertices_cff.ntracks
ntracks
Definition: vertices_cff.py:34
TransientTrack.h
edm::FileInPath
Definition: FileInPath.h:64
MakerMacros.h
testProducerWithPsetDescEmpty_cfi.x1
x1
Definition: testProducerWithPsetDescEmpty_cfi.py:33
edm::RefVector::end
const_iterator end() const
Termination of iteration.
Definition: RefVector.h:228
VtxSmearedMatchPbPBoost_cff.Beta
Beta
Definition: VtxSmearedMatchPbPBoost_cff.py:9
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
Provenance.h
JPTJetCollection.h
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
edm::ParameterSet
Definition: ParameterSet.h:47
Event.h
deltaR.h
LorentzVector.h
HE
Definition: GlobalHaloAlgo.cc:15
M_PI
#define M_PI
Definition: BXVectorInputProducer.cc:49
EE
Definition: GlobalHaloAlgo.cc:15
std
Definition: JetResolutionObject.h:76
metsig::jet
Definition: SignAlgoResolutions.h:47
edm::RefVectorIterator
Definition: EDProductfwd.h:33
reco::details::loadTMVAWeights
TMVA::IMethod * loadTMVAWeights(TMVA::Reader *reader, const std::string &method, const std::string &weightFile, bool verbose=false)
Definition: TMVAZipReader.cc:52
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
edm::RefToBase< reco::Jet >
TMVAZipReader.h
PileupJPTJetIdAlgo.h
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
ParameterSet.h
edm::RefVector::size
size_type size() const
Size of the RefVector.
Definition: RefVector.h:102
Vector3D.h
edm::FileInPath::fullPath
std::string fullPath() const
Definition: FileInPath.cc:163
JPTJet.h
cms
Namespace of DDCMS conversion namespace.
Definition: ProducerAnalyzer.cc:21
Generated for CMSSW Reference Manual by doxygen 1.8.16