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