CMS 3D CMS Logo

List of all members | Public Member Functions | Private Attributes
cms::PileupJPTJetIdAlgo Class Reference

#include <PileupJPTJetIdAlgo.h>

Public Member Functions

void bookMVAReader ()
 
float fillJPTBlock (const reco::JPTJet *jet)
 
 PileupJPTJetIdAlgo (const edm::ParameterSet &fParameters)
 
virtual ~PileupJPTJetIdAlgo ()
 

Private Attributes

float Beta
 
float dAxis1c
 
float dAxis1t
 
float dAxis2c
 
float dAxis2t
 
float EtaJ
 
float MultCalo
 
float MultTr
 
float Nvtx
 
float PtJ
 
TMVA::Reader * reader_
 
TMVA::Reader * readerF_
 
std::string tmvaMethod_
 
std::string tmvaWeights_
 
std::string tmvaWeightsF_
 
int verbosity
 

Detailed Description

Definition at line 27 of file PileupJPTJetIdAlgo.h.

Constructor & Destructor Documentation

◆ PileupJPTJetIdAlgo()

cms::PileupJPTJetIdAlgo::PileupJPTJetIdAlgo ( const edm::ParameterSet fParameters)

Definition at line 52 of file PileupJPTJetIdAlgo.cc.

References edm::FileInPath::fullPath(), edm::ParameterSet::getParameter(), AlCaHLTBitMon_QueryRunRegistry::string, and verbosity.

52  {
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  }
std::string fullPath() const
Definition: FileInPath.cc:161

◆ ~PileupJPTJetIdAlgo()

cms::PileupJPTJetIdAlgo::~PileupJPTJetIdAlgo ( )
virtual

Definition at line 59 of file PileupJPTJetIdAlgo.cc.

59  {
60  // std::cout<<" JetPlusTrack destructor "<<std::endl;
61  }

Member Function Documentation

◆ bookMVAReader()

void cms::PileupJPTJetIdAlgo::bookMVAReader ( )

Definition at line 63 of file PileupJPTJetIdAlgo.cc.

References VtxSmearedMatchPbPBoost_cff::Beta, gather_cfg::cout, reco::details::loadTMVAWeights(), and verbosity.

Referenced by PileupJPTJetIdProducer::PileupJPTJetIdProducer().

63  {
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 
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 
100 
101  // std::cout<<" Method booked F "<<std::endl;
102  }
TMVA::IMethod * loadTMVAWeights(TMVA::Reader *reader, const std::string &method, const std::string &weightFile, bool verbose=false)

◆ fillJPTBlock()

float cms::PileupJPTJetIdAlgo::fillJPTBlock ( const reco::JPTJet jet)

Definition at line 104 of file PileupJPTJetIdAlgo.cc.

References funct::abs(), edm::RefVector< C, T, F >::begin(), VtxSmearedMatchPbPBoost_cff::Beta, runTheMatrix::const, gather_cfg::cout, EE, edm::RefVector< C, T, F >::end(), HE, ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, metsig::jet, M_PI, vertices_cff::ntracks, edm::RefVector< C, T, F >::size(), mathSSE::sqrt(), and verbosity.

Referenced by PileupJPTJetIdProducer::produce().

104  {
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  const 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
Jets made from CaloTowers.
Definition: CaloJet.h:27
T sqrt(T t)
Definition: SSEVec.h:19
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
#define M_PI
size_type size() const
Size of the RefVector.
Definition: RefVector.h:102
const_iterator end() const
Termination of iteration.
Definition: RefVector.h:228
const_iterator begin() const
Initialize an iterator over the RefVector.
Definition: RefVector.h:223

Member Data Documentation

◆ Beta

float cms::PileupJPTJetIdAlgo::Beta
private

Definition at line 41 of file PileupJPTJetIdAlgo.h.

◆ dAxis1c

float cms::PileupJPTJetIdAlgo::dAxis1c
private

Definition at line 41 of file PileupJPTJetIdAlgo.h.

◆ dAxis1t

float cms::PileupJPTJetIdAlgo::dAxis1t
private

Definition at line 41 of file PileupJPTJetIdAlgo.h.

◆ dAxis2c

float cms::PileupJPTJetIdAlgo::dAxis2c
private

Definition at line 41 of file PileupJPTJetIdAlgo.h.

◆ dAxis2t

float cms::PileupJPTJetIdAlgo::dAxis2t
private

Definition at line 41 of file PileupJPTJetIdAlgo.h.

◆ EtaJ

float cms::PileupJPTJetIdAlgo::EtaJ
private

Definition at line 41 of file PileupJPTJetIdAlgo.h.

◆ MultCalo

float cms::PileupJPTJetIdAlgo::MultCalo
private

Definition at line 41 of file PileupJPTJetIdAlgo.h.

◆ MultTr

float cms::PileupJPTJetIdAlgo::MultTr
private

Definition at line 41 of file PileupJPTJetIdAlgo.h.

◆ Nvtx

float cms::PileupJPTJetIdAlgo::Nvtx
private

Definition at line 41 of file PileupJPTJetIdAlgo.h.

◆ PtJ

float cms::PileupJPTJetIdAlgo::PtJ
private

Definition at line 41 of file PileupJPTJetIdAlgo.h.

◆ reader_

TMVA::Reader* cms::PileupJPTJetIdAlgo::reader_
private

Definition at line 42 of file PileupJPTJetIdAlgo.h.

◆ readerF_

TMVA::Reader* cms::PileupJPTJetIdAlgo::readerF_
private

Definition at line 43 of file PileupJPTJetIdAlgo.h.

◆ tmvaMethod_

std::string cms::PileupJPTJetIdAlgo::tmvaMethod_
private

Definition at line 44 of file PileupJPTJetIdAlgo.h.

◆ tmvaWeights_

std::string cms::PileupJPTJetIdAlgo::tmvaWeights_
private

Definition at line 44 of file PileupJPTJetIdAlgo.h.

◆ tmvaWeightsF_

std::string cms::PileupJPTJetIdAlgo::tmvaWeightsF_
private

Definition at line 44 of file PileupJPTJetIdAlgo.h.

◆ verbosity

int cms::PileupJPTJetIdAlgo::verbosity
private

Definition at line 38 of file PileupJPTJetIdAlgo.h.