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 28 of file PileupJPTJetIdAlgo.h.

Constructor & Destructor Documentation

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

Definition at line 54 of file PileupJPTJetIdAlgo.cc.

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

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  }
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
edm::FileInPath::fullPath
std::string fullPath() const
Definition: FileInPath.cc:184
cms::PileupJPTJetIdAlgo::~PileupJPTJetIdAlgo ( )
virtual

Definition at line 64 of file PileupJPTJetIdAlgo.cc.

65  {
66  // std::cout<<" JetPlusTrack destructor "<<std::endl;
67  }

Member Function Documentation

void cms::PileupJPTJetIdAlgo::bookMVAReader ( )

Definition at line 69 of file PileupJPTJetIdAlgo.cc.

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

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  }
reco::details::loadTMVAWeights
TMVA::IMethod * loadTMVAWeights(TMVA::Reader *reader, const std::string &method, const std::string &weightFile, bool verbose=false)
Definition: TMVAZipReader.cc:57
gather_cfg.cout
cout
Definition: gather_cfg.py:145
float cms::PileupJPTJetIdAlgo::fillJPTBlock ( const reco::JPTJet *  jet)

Definition at line 112 of file PileupJPTJetIdAlgo.cc.

References funct::abs(), edm::RefVector< C, T, F >::begin(), gather_cfg::cout, EE, edm::RefVector< C, T, F >::end(), reco::JPTJet::getCaloJetRef(), HE, M_PI, edm::RefVector< C, T, F >::size(), mathSSE::sqrt(), and HIPAlignmentAlgorithm_cfi::verbosity.

Referenced by PileupJPTJetIdProducer::produce().

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 EE=0.;
134  double HE=0.;
135  double EELong=0.;
136  double EEShort=0.;
137 
138  std::vector <CaloTowerPtr> calotwrs = (*rawcalojet).getCaloConstituents();
139 
140  if (verbosity > 0) { std::cout<<"======= CaloTowerPtr DONE == " << std::endl;}
141 
142  for(std::vector <CaloTowerPtr>::const_iterator icalot = calotwrs.begin();
143  icalot!= calotwrs.end(); icalot++) {
144  ncalotowers++;
145 
146  double deta=(*jet).eta()-(*icalot)->eta();
147  double dphi=(*jet).phi()-(*icalot)->phi();
148 
149  if(dphi > M_PI ) dphi = dphi-2.*M_PI;
150  if(dphi < -1.*M_PI ) dphi = dphi+2.*M_PI;
151 
152  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;
153 
154 
155  if(abs((*icalot)->ieta())<30) EE = EE + (*icalot)->emEnergy();
156  if(abs((*icalot)->ieta())<30) HE = HE + (*icalot)->hadEnergy();
157  if(abs((*icalot)->ieta())>29) EELong = EELong + (*icalot)->emEnergy();
158  if(abs((*icalot)->ieta())>29) EEShort = EEShort + (*icalot)->hadEnergy();
159 
160  sumpt = sumpt + (*icalot)->pt();
161  dphi2 = dphi2 + dphi*dphi*(*icalot)->pt();
162  deta2 = deta2 + deta*deta*(*icalot)->pt();
163  dphi1 = dphi1 + dphi*(*icalot)->pt();
164  deta1 = deta1 + deta*(*icalot)->pt();
165  dphideta = dphideta + dphi*deta*(*icalot)->pt();
166 
167  } // calojet constituents
168 
169 
170  if( sumpt > 0.) {
171  deta1 = deta1/sumpt;
172  dphi1 = dphi1/sumpt;
173  deta2 = deta2/sumpt;
174  dphi2 = dphi2/sumpt;
175  dphideta = dphideta/sumpt;
176  }
177 
178  // W.r.t. principal axis
179 
180  double detavar = deta2-deta1*deta1;
181  double dphivar = dphi2-dphi1*dphi1;
182  double dphidetacov = dphideta - deta1*dphi1;
183 
184  double det = (detavar-dphivar)*(detavar-dphivar)+4*dphidetacov*dphidetacov;
185  det = sqrt(det);
186  double x1 = (detavar+dphivar+det)/2.;
187  double x2 = (detavar+dphivar-det)/2.;
188 
189 
190 
191 if (verbosity > 0)
192 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;
193 
194 
195 // For jets with |eta|<2 take also tracks shape
196  int ntracks=0.;
197  double sumpttr=0.;
198  double dphitr2=0.;
199  double detatr2=0.;
200  double dphitr1=0.;
201  double detatr1=0.;
202  double dphidetatr=0.;
203 
204  const reco::TrackRefVector pioninin = (*jet).getPionsInVertexInCalo();
205 
206  for(reco::TrackRefVector::const_iterator it = pioninin.begin(); it != pioninin.end(); it++) {
207  if ((*it)->pt() > 0.5 && ((*it)->ptError()/(*it)->pt()) < 0.05 )
208  {
209  ntracks++;
210  sumpttr = sumpttr + (*it)->pt();
211  double deta=(*jet).eta()-(*it)->eta();
212  double dphi=(*jet).phi()-(*it)->phi();
213 
214  if(dphi > M_PI ) dphi = dphi-2.*M_PI;
215  if(dphi < -1.*M_PI ) dphi = dphi+2.*M_PI;
216 
217  dphitr2 = dphitr2 + dphi*dphi*(*it)->pt();
218  detatr2 = detatr2 + deta*deta*(*it)->pt();
219  dphitr1 = dphitr1 + dphi*(*it)->pt();
220  detatr1 = detatr1 + deta*(*it)->pt();
221  dphidetatr = dphidetatr + dphi*deta*(*it)->pt();
222  if(verbosity>0) std::cout<<" Tracks-in-in "<<(*it)->pt()<<" "<<(*it)->eta()<<" "<<(*it)->phi()<<" in jet "<<(*jet).eta()<<" "<<
223  (*jet).phi()<<" jet pt "<<(*jet).pt()<<std::endl;
224  }
225  }// pioninin
226 
227 
228  const reco::TrackRefVector pioninout = (*jet).getPionsInVertexOutCalo();
229 
230  for(reco::TrackRefVector::const_iterator it = pioninout.begin(); it != pioninout.end(); it++) {
231  if ((*it)->pt() > 0.5 && ((*it)->ptError()/(*it)->pt()) < 0.05 )
232  {
233  ntracks++;
234  sumpttr = sumpttr + (*it)->pt();
235  double deta=(*jet).eta()-(*it)->eta();
236  double dphi=(*jet).phi()-(*it)->phi();
237 
238  if(dphi > M_PI ) dphi = dphi-2.*M_PI;
239  if(dphi < -1.*M_PI ) dphi = dphi+2.*M_PI;
240 
241  dphitr2 = dphitr2 + dphi*dphi*(*it)->pt();
242  detatr2 = detatr2 + deta*deta*(*it)->pt();
243  dphitr1 = dphitr1 + dphi*(*it)->pt();
244  detatr1 = detatr1 + deta*(*it)->pt();
245  dphidetatr = dphidetatr + dphi*deta*(*it)->pt();
246  if(verbosity>0) std::cout<<" Tracks-in-in "<<(*it)->pt()<<" "<<(*it)->eta()<<" "<<(*it)->phi()<<" in jet "<<(*jet).eta()<<" "<<
247  (*jet).phi()<<" jet pt "<<(*jet).pt()<<std::endl;
248  }
249  }// pioninout
250 
251  if(verbosity>0) std::cout<<" Number of tracks in-in and in-out "<<pioninin.size()<<" "<<pioninout.size()<<std::endl;
252 
253 
254  if( sumpttr > 0.) {
255  detatr1 = detatr1/sumpttr;
256  dphitr1 = dphitr1/sumpttr;
257  detatr2 = detatr2/sumpttr;
258  dphitr2 = dphitr2/sumpttr;
259  dphidetatr = dphidetatr/sumpttr;
260  }
261 
262  // W.r.t. principal axis
263 
264  double detavart = detatr2-detatr1*detatr1;
265  double dphivart = dphitr2-dphitr1*dphitr1;
266  double dphidetacovt = dphidetatr - detatr1*dphitr1;
267 
268  double dettr = (detavart-dphivart)*(detavart-dphivart)+4*dphidetacovt*dphidetacovt;
269  dettr = sqrt(dettr);
270  double x1tr = (detavart+dphivart+dettr)/2.;
271  double x2tr = (detavart+dphivart-dettr)/2.;
272 
273  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;
274 
275 // Multivariate analisis
276  PtJ = (*jet).pt();
277  EtaJ = (*jet).eta();
278  Beta = (*jet).getSpecific().Zch;
279  dAxis2c = x2;
280  dAxis1c = x1;
281  dAxis2t = x2tr;
282  dAxis1t = x1tr;
283  MultCalo = ncalotowers;
284  MultTr = ntracks;
285 
286  float mva_ = 1.;
287  if(fabs(EtaJ)<2.6) {
288  mva_ = reader_->EvaluateMVA( "BDTG method" );
289  // std::cout<<" MVA analysis "<<mva_<<std::endl;
290  } else {
291  mva_ = readerF_->EvaluateMVA( "BDTG method" );
292  // std::cout<<" MVA analysis Forward "<<mva_<<std::endl;
293  }
294  if (verbosity > 0) std::cout<<"======= Computed MVA = " <<
295  mva_ <<std::endl;
296  return mva_;
297 } // fillJPTBlock
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:253
edm::RefVector::begin
const_iterator begin() const
Initialize an iterator over the RefVector.
Definition: RefVector.h:248
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
M_PI
#define M_PI
Definition: BXVectorInputProducer.cc:51
edm::RefVector::size
size_type size() const
Size of the RefVector.
Definition: RefVector.h:107
gather_cfg.cout
cout
Definition: gather_cfg.py:145

Member Data Documentation

float cms::PileupJPTJetIdAlgo::Beta
private

Definition at line 44 of file PileupJPTJetIdAlgo.h.

float cms::PileupJPTJetIdAlgo::dAxis1c
private

Definition at line 44 of file PileupJPTJetIdAlgo.h.

float cms::PileupJPTJetIdAlgo::dAxis1t
private

Definition at line 44 of file PileupJPTJetIdAlgo.h.

float cms::PileupJPTJetIdAlgo::dAxis2c
private

Definition at line 44 of file PileupJPTJetIdAlgo.h.

float cms::PileupJPTJetIdAlgo::dAxis2t
private

Definition at line 44 of file PileupJPTJetIdAlgo.h.

float cms::PileupJPTJetIdAlgo::EtaJ
private

Definition at line 44 of file PileupJPTJetIdAlgo.h.

float cms::PileupJPTJetIdAlgo::MultCalo
private

Definition at line 44 of file PileupJPTJetIdAlgo.h.

float cms::PileupJPTJetIdAlgo::MultTr
private

Definition at line 44 of file PileupJPTJetIdAlgo.h.

float cms::PileupJPTJetIdAlgo::Nvtx
private

Definition at line 44 of file PileupJPTJetIdAlgo.h.

float cms::PileupJPTJetIdAlgo::PtJ
private

Definition at line 44 of file PileupJPTJetIdAlgo.h.

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

Definition at line 45 of file PileupJPTJetIdAlgo.h.

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

Definition at line 46 of file PileupJPTJetIdAlgo.h.

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

Definition at line 47 of file PileupJPTJetIdAlgo.h.

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

Definition at line 47 of file PileupJPTJetIdAlgo.h.

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

Definition at line 47 of file PileupJPTJetIdAlgo.h.

int cms::PileupJPTJetIdAlgo::verbosity
private

Definition at line 41 of file PileupJPTJetIdAlgo.h.

Generated for CMSSW Reference Manual by doxygen 1.8.11