| 1 |
// -*- C++ -*-
|
| 2 |
//
|
| 3 |
// Package: PATJetIDAnalyzer
|
| 4 |
// Class: PATJetIDAnalyzer
|
| 5 |
//
|
| 6 |
/**\class PATJetIDAnalyzer PATJetIDAnalyzer.cc Demo/PATJetIDAnalyzer/src/PATJetIDAnalyzer.cc
|
| 7 |
|
| 8 |
Description: <one line class summary>
|
| 9 |
|
| 10 |
Implementation:
|
| 11 |
<Notes on implementation>
|
| 12 |
*/
|
| 13 |
//
|
| 14 |
// Original Author: Christian Autermann
|
| 15 |
// Created: Mon Feb 25 11:33:02 CET 2008
|
| 16 |
// $Id: PATJetIDAnalyzer.cc,v 1.2 2008/02/25 17:52:29 auterman Exp $
|
| 17 |
//
|
| 18 |
//
|
| 19 |
|
| 20 |
//system includes
|
| 21 |
#include <fstream>
|
| 22 |
#include <iostream>
|
| 23 |
//framework includes
|
| 24 |
#include "FWCore/Utilities/interface/EDMException.h"
|
| 25 |
#include "FWCore/ServiceRegistry/interface/Service.h"
|
| 26 |
#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
|
| 27 |
//root includes
|
| 28 |
#include "TH1F.h"
|
| 29 |
#include "TH2F.h"
|
| 30 |
#include "TString.h"
|
| 31 |
#include "TRandom.h"
|
| 32 |
//User include files
|
| 33 |
#include "DataFormats/Math/interface/deltaR.h"
|
| 34 |
//AOD
|
| 35 |
#include "DataFormats/JetReco/interface/GenJet.h"
|
| 36 |
#include "DataFormats/JetReco/interface/GenJetfwd.h"
|
| 37 |
#include "DataFormats/JetReco/interface/CaloJetCollection.h"
|
| 38 |
#include "DataFormats/METReco/interface/CaloMETCollection.h"
|
| 39 |
#include "DataFormats/METReco/interface/GenMETCollection.h"
|
| 40 |
#include "DataFormats/HepMCCandidate/interface/GenParticleCandidate.h"
|
| 41 |
//PAT
|
| 42 |
#include "DataFormats/PatCandidates/interface/Jet.h"
|
| 43 |
#include "DataFormats/PatCandidates/interface/MET.h"
|
| 44 |
|
| 45 |
#include "PhysicsTools/Utilities/interface/EtComparator.h"
|
| 46 |
#include "Demo/PATJetIDAnalyzer/interface/NameScheme.h"
|
| 47 |
#include "Demo/PATJetIDAnalyzer/interface/PATJetIDAnalyzer.h"
|
| 48 |
#include "TMatrixT.h"
|
| 49 |
|
| 50 |
//EM-cells
|
| 51 |
#include "DataFormats/EcalDetId/interface/EBDetId.h"
|
| 52 |
#include "Geometry/Records/interface/IdealGeometryRecord.h"
|
| 53 |
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
|
| 54 |
#include "FWCore/Framework/interface/ESHandle.h"
|
| 55 |
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
|
| 56 |
#include "FWCore/Framework/interface/EventSetup.h"
|
| 57 |
|
| 58 |
|
| 59 |
const unsigned PATJetIDAnalyzer::N_Fourier_Bins_2D;
|
| 60 |
const unsigned PATJetIDAnalyzer::N_Fourier_Bins_1D;
|
| 61 |
const unsigned PATJetIDAnalyzer::Nf_Fourier_Bins_1D;
|
| 62 |
//
|
| 63 |
// constructors and destructor
|
| 64 |
//
|
| 65 |
PATJetIDAnalyzer::PATJetIDAnalyzer(const edm::ParameterSet& iConfig) :
|
| 66 |
_patJet ( iConfig.getParameter<edm::InputTag>( "patJet" ) ),
|
| 67 |
_patMet ( iConfig.getParameter<edm::InputTag>( "patMet" ) ),
|
| 68 |
_ebrechits( iConfig.getParameter<edm::InputTag>("EBRecHits") ),
|
| 69 |
_hist ( iConfig.getParameter<std::string>( "hist" ) ),
|
| 70 |
_jetminpt ( iConfig.getParameter<double>( "MinJetPt" ) ),
|
| 71 |
_jetmaxeta(iConfig.getParameter<double>( "MaxJetEta" ) ),
|
| 72 |
_simulate_noise(iConfig.getParameter<bool>( "SimulateNoise" ) ),
|
| 73 |
_NoiseMean(iConfig.getParameter<double>( "NoiseMean" ) ),
|
| 74 |
_NoiseSigma(iConfig.getParameter<double>( "NoiseSigma" ) ),
|
| 75 |
_NoiseThreshold(iConfig.getParameter<double>( "NoiseThreshold" ) ),
|
| 76 |
_DoNormalization(iConfig.getParameter<bool>( "Normalize" ) ),
|
| 77 |
_uniqueplotid(0)
|
| 78 |
{
|
| 79 |
//now do what ever initialization is needed
|
| 80 |
random = new TRandom(123455);
|
| 81 |
}
|
| 82 |
|
| 83 |
|
| 84 |
PATJetIDAnalyzer::~PATJetIDAnalyzer()
|
| 85 |
{
|
| 86 |
|
| 87 |
// do anything here that needs to be done at desctruction time
|
| 88 |
// (e.g. close files, deallocate resources etc.)
|
| 89 |
|
| 90 |
}
|
| 91 |
|
| 92 |
|
| 93 |
//
|
| 94 |
// member functions
|
| 95 |
//
|
| 96 |
|
| 97 |
double PATJetIDAnalyzer::norm(int n)
|
| 98 |
{
|
| 99 |
return ( n==0 ? 1 : sqrt(2.) );
|
| 100 |
}
|
| 101 |
|
| 102 |
TH2F * PATJetIDAnalyzer::FDCT(TH2F & k)
|
| 103 |
{
|
| 104 |
TH2F* result= new TH2F(k);
|
| 105 |
unsigned int Nx=k.GetNbinsX();
|
| 106 |
unsigned int Ny=k.GetNbinsY();
|
| 107 |
|
| 108 |
for (unsigned int u=0; u<=Nx; ++u){
|
| 109 |
for (unsigned int v=0; v<=Ny; ++v){
|
| 110 |
result->SetBinContent(u,v,0.0);
|
| 111 |
for (unsigned int i=0; i<Nx; ++i){
|
| 112 |
for (unsigned int j=0; j<Ny; ++j){
|
| 113 |
result->SetBinContent(u,v, result->GetBinContent(u,v)+
|
| 114 |
norm(u)*norm(v)/(sqrt((double)Nx)*sqrt((double)Ny))*
|
| 115 |
(k.GetBinContent(i,j)-0.)*
|
| 116 |
cos( ( (double)i+0.5)*3.14159*(double)u/( (double)Nx) )*
|
| 117 |
cos( ( (double)j+0.5)*3.14159*(double)v/( (double)Ny) )
|
| 118 |
);
|
| 119 |
//cout << "x="<<u<<", y="<<v<<", f="<<result->GetBinContent(u,v)<<endl;
|
| 120 |
}
|
| 121 |
}
|
| 122 |
}
|
| 123 |
}
|
| 124 |
return result;
|
| 125 |
}
|
| 126 |
|
| 127 |
void PATJetIDAnalyzer::FakeNoise(TH2F & hk, double * k)
|
| 128 |
{
|
| 129 |
unsigned int Nx=hk.GetNbinsX();
|
| 130 |
unsigned int Ny=hk.GetNbinsY();
|
| 131 |
unsigned int r;
|
| 132 |
double noise, absnoise=0., abs=0.;
|
| 133 |
double norm = hk.Integral();
|
| 134 |
for (unsigned i=0; i<N_Fourier_Bins_1D; ++i)
|
| 135 |
abs+=k[i];
|
| 136 |
for (unsigned int u=1; u<=Nx; ++u){
|
| 137 |
for (unsigned int v=1; v<=Ny; ++v){
|
| 138 |
noise = random->Gaus(_NoiseMean,_NoiseSigma);
|
| 139 |
if (noise+hk.GetBinContent(u,v)>_NoiseThreshold) {
|
| 140 |
hk.SetBinContent(u,v, noise+hk.GetBinContent(u,v) );
|
| 141 |
double phi = ((double)(u-1)-(double)Nx)/2/(double)Nx;
|
| 142 |
double eta = ((double)(v-1)-(double)Ny)/2/(double)Ny;
|
| 143 |
r = (int)(sqrt( (phi*phi+eta*eta)/2. )*(double)(N_Fourier_Bins_1D-1)/0.5);
|
| 144 |
//cout <<"u="<<u<<", v="<<v<<", phi="<<phi<<", eta="<<eta<<", r="<<r<<", E(r)="<<k[r]
|
| 145 |
// <<", E(u,v)="<<hk.GetBinContent(u,v)<<", noise="<<noise<<endl;
|
| 146 |
if (r>0 && r<N_Fourier_Bins_1D){ k[r] += noise; absnoise+=noise; }
|
| 147 |
}
|
| 148 |
}
|
| 149 |
}
|
| 150 |
|
| 151 |
//normalize such, that the energy of the jet ("brightness") is unchanged by noise
|
| 152 |
if (_DoNormalization){
|
| 153 |
hk.Scale( norm/hk.Integral() );
|
| 154 |
for (unsigned i=0; i<N_Fourier_Bins_1D; ++i)
|
| 155 |
k[i]=k[i]*abs/(abs+absnoise);
|
| 156 |
}
|
| 157 |
}
|
| 158 |
|
| 159 |
void PATJetIDAnalyzer::FakeNoise2x2(TH2F & hk, double * k)
|
| 160 |
{
|
| 161 |
unsigned int Nx=hk.GetNbinsX();
|
| 162 |
unsigned int Ny=hk.GetNbinsY();
|
| 163 |
unsigned int r;
|
| 164 |
double noise, absnoise=0., abs=0.;
|
| 165 |
double norm = hk.Integral();
|
| 166 |
for (unsigned i=0; i<N_Fourier_Bins_1D; ++i)
|
| 167 |
abs+=k[i];
|
| 168 |
unsigned int u=(int)(random->Rndm()*(double)hk.GetNbinsX())+1;
|
| 169 |
unsigned int v=(int)(random->Rndm()*(double)hk.GetNbinsY())+1;
|
| 170 |
noise = 15.;
|
| 171 |
cout << u << ", " << v << ": old="<<hk.GetBinContent(u,v);
|
| 172 |
hk.SetBinContent(u,v, noise+hk.GetBinContent(u,v) );
|
| 173 |
cout << "; new="<<hk.GetBinContent(u,v)<<endl;;
|
| 174 |
double phi = ((double)(u-1)-(double)Nx)/2/(double)Nx;
|
| 175 |
double eta = ((double)(v-1)-(double)Ny)/2/(double)Ny;
|
| 176 |
r = (int)(sqrt( (phi*phi+eta*eta)/2. )*(double)(N_Fourier_Bins_1D-1)/0.5);
|
| 177 |
//cout <<"u="<<u<<", v="<<v<<", phi="<<phi<<", eta="<<eta<<", r="<<r<<", E(r)="<<k[r]
|
| 178 |
// <<", E(u,v)="<<hk.GetBinContent(u,v)<<", noise="<<noise<<endl;
|
| 179 |
if (r>0 && r<N_Fourier_Bins_1D){ k[r] += noise; absnoise+=noise; }
|
| 180 |
|
| 181 |
//normalize such, that the energy of the jet ("brightness") is unchanged by noise
|
| 182 |
if (_DoNormalization){
|
| 183 |
hk.Scale( norm/hk.Integral() );
|
| 184 |
for (unsigned i=0; i<N_Fourier_Bins_1D; ++i)
|
| 185 |
k[i]=k[i]*abs/(abs+absnoise);
|
| 186 |
}
|
| 187 |
}
|
| 188 |
|
| 189 |
double PATJetIDAnalyzer::dphi(double phi1, double phi2)
|
| 190 |
{
|
| 191 |
double result = phi1 - phi2;
|
| 192 |
if (result> 3.14159) result-=3.14159;
|
| 193 |
if (result<-3.14159) result+=3.14159;
|
| 194 |
return result;
|
| 195 |
}
|
| 196 |
|
| 197 |
|
| 198 |
void PATJetIDAnalyzer::FourierTransformation( const unsigned int i,
|
| 199 |
const pat::Jet& jet,
|
| 200 |
const EBRecHitCollection& EBRecHit)
|
| 201 |
{
|
| 202 |
///Studies
|
| 203 |
double k[N_Fourier_Bins_1D], f[Nf_Fourier_Bins_1D];
|
| 204 |
for (unsigned b=0; b<N_Fourier_Bins_1D; ++b){
|
| 205 |
k[b]=0.;
|
| 206 |
//k[b] = cos(2.*3.1415*(double)b/(double)N_Fourier_Bins) +
|
| 207 |
// cos(0.2*2.*3.1415*(double)b/(double)N_Fourier_Bins) +
|
| 208 |
// cos(10*2.*3.1415*(double)b/(double)N_Fourier_Bins);
|
| 209 |
|
| 210 |
}
|
| 211 |
for (unsigned b=0; b<Nf_Fourier_Bins_1D; ++b){
|
| 212 |
f[b]=0.;
|
| 213 |
}
|
| 214 |
char * name = new char[64];
|
| 215 |
sprintf(name,"h_%d",++_uniqueplotid);
|
| 216 |
TH2F * k2d = new TH2F(name,"",16,-8,7,16,-8,7);
|
| 217 |
TH2F * f2d;
|
| 218 |
|
| 219 |
int jtow=0; //, icell=0;
|
| 220 |
std::vector <CaloTowerRef> jetTowers = jet.getConstituents();
|
| 221 |
|
| 222 |
for(std::vector<CaloTowerRef>::const_iterator tow = jetTowers.begin();
|
| 223 |
tow != jetTowers.end(); ++tow, ++jtow){
|
| 224 |
|
| 225 |
//2D-energy picture of this jet
|
| 226 |
k2d->Fill( dphi((*tow)->phi(),jet.phi())/0.0873,
|
| 227 |
((*tow)->eta()-jet.eta())/0.087,
|
| 228 |
(*tow)->energy() );
|
| 229 |
|
| 230 |
//1D-energy picture
|
| 231 |
double dR = deltaR(jet,*(*tow));
|
| 232 |
if (fabs(dR)>0.5) dR=0.5;
|
| 233 |
//_ft_k[i]->Fill( dR, (*tow)->energy() );
|
| 234 |
k[ (int)(((double)N_Fourier_Bins_1D-1.)*dR/0.5) ] += (*tow)->energy();
|
| 235 |
|
| 236 |
//no EM cells in my current test data....
|
| 237 |
/*
|
| 238 |
double eem=0.;
|
| 239 |
for (size_t it=0; it<(*tow)->constituentsSize(); ++it) {
|
| 240 |
const DetId detid = (*tow)->constituent(it);
|
| 241 |
EcalRecHitCollection::const_iterator myRecHit = EBRecHit.find(detid);
|
| 242 |
if(myRecHit != EBRecHit.end()) {
|
| 243 |
eem += myRecHit->energy();
|
| 244 |
EBDetId det = myRecHit->id();
|
| 245 |
|
| 246 |
const CaloCellGeometry* cell=EBgeom->getGeometry( myRecHit->id() );
|
| 247 |
/
|
| 248 |
etowet [icell] = myRecHit->energy()*sin( cell->getPosition().theta());
|
| 249 |
etoweta[icell] = cell->getPosition().eta();
|
| 250 |
etowphi[icell] = cell->getPosition().phi();
|
| 251 |
etowe [icell] = myRecHit->energy();
|
| 252 |
etowid_phi[icell] = det.iphi();
|
| 253 |
etowid_eta[icell] = det.ieta();
|
| 254 |
etowid [icell] = myRecHit->id().rawId();
|
| 255 |
etownum[icell] = icell;
|
| 256 |
//
|
| 257 |
_ft_energy[i]->Fill( deltaPhi( (double)cell->getPosition().phi(),jet.phi()), cell->getPosition().eta()-jet.eta(), myRecHit->energy() );
|
| 258 |
++icell;
|
| 259 |
}
|
| 260 |
}
|
| 261 |
*/
|
| 262 |
}
|
| 263 |
|
| 264 |
|
| 265 |
//Add fake-noise to the jet:
|
| 266 |
if(_simulate_noise)
|
| 267 |
FakeNoise( *k2d, k);
|
| 268 |
//FakeNoise2x2( *k2d, k);
|
| 269 |
|
| 270 |
|
| 271 |
//1-D fourier transformation
|
| 272 |
for (unsigned b=0; b<N_Fourier_Bins_1D; ++b) {
|
| 273 |
_ft_k[i]->SetBinContent(b, k[b] + _ft_k[i]->GetBinContent(b) );
|
| 274 |
//before fourier-trafo substract average jet spectrum:
|
| 275 |
k[b] -= 167.53*exp( -11.941*0.5*(double)b/(double)N_Fourier_Bins_1D );
|
| 276 |
_ft_ksubavg[i]->SetBinContent(b, k[b] + _ft_ksubavg[i]->GetBinContent(b) );
|
| 277 |
}
|
| 278 |
for (unsigned b=0; b<=Nf_Fourier_Bins_1D; ++b) {
|
| 279 |
double fourier=0.;
|
| 280 |
for (unsigned j=0; j<N_Fourier_Bins_1D; ++j) {
|
| 281 |
fourier += 2./0.5*k[j]*
|
| 282 |
cos( ((double)j)*3.14159*(double)b/
|
| 283 |
(double)N_Fourier_Bins_1D );
|
| 284 |
}
|
| 285 |
//_ft_f[i]->SetBinContent(b, fourier);
|
| 286 |
_ft_f[i]->SetBinContent(b+1, fourier + _ft_f[i]->GetBinContent(b+1) );
|
| 287 |
}
|
| 288 |
|
| 289 |
|
| 290 |
//2D-fourier transformation & calculation of observables
|
| 291 |
TMatrixT<double> matrix16(16,16);
|
| 292 |
TMatrixT<double> matrix8(8,8);
|
| 293 |
for (int binx=0; binx<=k2d->GetNbinsX();++binx) {
|
| 294 |
for (int biny=0; biny<=k2d->GetNbinsY();++biny) {
|
| 295 |
_ft_energy[i]->SetBinContent(binx,biny,_ft_energy[i]->GetBinContent(binx,biny)+
|
| 296 |
k2d->GetBinContent(binx,biny) );
|
| 297 |
}}
|
| 298 |
|
| 299 |
std::vector<double> db;
|
| 300 |
double sum=0.;
|
| 301 |
f2d = FDCT( *k2d);
|
| 302 |
for (int binx=0; binx<=f2d->GetNbinsX();++binx) {
|
| 303 |
for (int biny=0; biny<=f2d->GetNbinsY();++biny) {
|
| 304 |
double bincont = f2d->GetBinContent(binx,biny);
|
| 305 |
_ft_frequency[i]->SetBinContent(binx,biny, _ft_frequency[i]->GetBinContent(binx,biny)+
|
| 306 |
bincont );
|
| 307 |
if (bincont!=0.) db.push_back( fabs(bincont) );
|
| 308 |
sum += fabs(bincont);
|
| 309 |
if (binx>0 && binx<16 &&biny>0 && biny<16)
|
| 310 |
matrix16[binx][biny] = f2d->GetBinContent(binx,biny);
|
| 311 |
if (binx>3 && binx<12 &&biny>3 && biny<12) {
|
| 312 |
matrix8[binx-4][biny-4] = f2d->GetBinContent(binx,biny);
|
| 313 |
//cout <<"(" <<binx << "::"<<biny<< ") = " <<k2d->GetBinContent(bin)<<endl;
|
| 314 |
}
|
| 315 |
}}
|
| 316 |
std::sort(db.begin(), db.end());
|
| 317 |
double ig=0., F10=0., Fs10=0.;
|
| 318 |
int n10=-1, n30=-1, n60=-1, n90=-1, n95=-1, n99=-1;
|
| 319 |
for (unsigned j=0; j<db.size(); ++j){
|
| 320 |
if ( ig>0.1*sum && n10==-1) n10 = j;
|
| 321 |
if ( ig>0.3*sum && n30==-1) n30 = j;
|
| 322 |
if ( ig>0.6*sum && n60==-1) n60 = j;
|
| 323 |
if ( ig>0.9*sum && n90==-1) n90 = j;
|
| 324 |
if ( ig>0.95*sum && n95==-1) n95 = j;
|
| 325 |
if ( ig>0.99*sum && n99==-1) n99 = j;
|
| 326 |
ig += fabs(db[j]);
|
| 327 |
if (j<10) Fs10 += db[j];
|
| 328 |
if (j>=db.size()-10 && db.size()>10) F10 += db[j];
|
| 329 |
}
|
| 330 |
double Ilow = f2d->Integral( 1, (int)(f2d->GetNbinsX()/3.),0,(int)(f2d->GetNbinsY()/3.) );
|
| 331 |
double Ihi = f2d->Integral( (int)(f2d->GetNbinsX()*2./3.),f2d->GetNbinsX(), (int)(f2d->GetNbinsY()*2./3.), f2d->GetNbinsY()-1 );
|
| 332 |
|
| 333 |
//fill observables from the 2D-fourier-transformation:
|
| 334 |
_fto_n99[i]->Fill( n99 );
|
| 335 |
_fto_n95[i]->Fill( n95 );
|
| 336 |
_fto_n90[i]->Fill( n90 );
|
| 337 |
_fto_n60[i]->Fill( n60 );
|
| 338 |
_fto_n30[i]->Fill( n30 );
|
| 339 |
_fto_n10[i]->Fill( n10 );
|
| 340 |
_fto_F10[i]->Fill( F10 );
|
| 341 |
_fto_Fs10[i]->Fill( Fs10 );
|
| 342 |
if (Ihi!=0.) {
|
| 343 |
_fto_LowFovHiF[i]->Fill( Ilow/Ihi );
|
| 344 |
_fto_LowFvsHiF[i]->Fill( (Ihi-Ilow)/(Ihi+Ilow) );
|
| 345 |
}
|
| 346 |
//_fto_det16[i]->Fill( matrix16.Determinant() );
|
| 347 |
//_fto_det8[i]->Fill( matrix8.Determinant() );
|
| 348 |
|
| 349 |
//clean-up
|
| 350 |
delete name;
|
| 351 |
delete f2d;
|
| 352 |
delete k2d;
|
| 353 |
}
|
| 354 |
|
| 355 |
// ------------ method called to for each event ------------
|
| 356 |
void
|
| 357 |
PATJetIDAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
|
| 358 |
{
|
| 359 |
using namespace edm;
|
| 360 |
using namespace reco;
|
| 361 |
using namespace std;
|
| 362 |
using namespace pat;
|
| 363 |
|
| 364 |
// edm::ESHandle<CaloGeometry> pG;
|
| 365 |
// iSetup.get<IdealGeometryRecord>().get(pG);
|
| 366 |
// const CaloGeometry cG = *pG;
|
| 367 |
// EBgeom=cG.getSubdetectorGeometry(DetId::Ecal,1);
|
| 368 |
|
| 369 |
//EB RecHits
|
| 370 |
const EBRecHitCollection *EBRecHit = 0;
|
| 371 |
edm::Handle<EBRecHitCollection> EcalRecHitEB;
|
| 372 |
// iEvent.getByLabel( _ebrechits, EcalRecHitEB);
|
| 373 |
// if( EcalRecHitEB.isValid() ){
|
| 374 |
// EBRecHit = EcalRecHitEB.product();
|
| 375 |
// }
|
| 376 |
|
| 377 |
//PAT Jets
|
| 378 |
typedef vector<pat::Jet>::const_iterator PatJetCIter;
|
| 379 |
edm::Handle<vector<pat::Jet> > PatJets;
|
| 380 |
iEvent.getByLabel( _patJet, PatJets );
|
| 381 |
GreaterByEt<pat::Jet> eTComparator_;
|
| 382 |
vector<pat::Jet> patjets = *PatJets;
|
| 383 |
std::sort(patjets.begin(), patjets.end(), eTComparator_);
|
| 384 |
PatJetCIter PatJetBegin = patjets.begin();
|
| 385 |
PatJetCIter PatJetEnd = patjets.end();
|
| 386 |
|
| 387 |
//Jets
|
| 388 |
unsigned multiplicity = 0;
|
| 389 |
double met=0.0, ht=0.0;
|
| 390 |
for (PatJetCIter it=PatJetBegin; it!=PatJetEnd; ++it) {
|
| 391 |
if (it->pt() > _jetminpt && fabs( it->eta() )<_jetmaxeta){
|
| 392 |
if (multiplicity<_njets ) {
|
| 393 |
_pt_jet[ multiplicity]->Fill( it->pt() );
|
| 394 |
_eta_jet[ multiplicity]->Fill( it->eta() );
|
| 395 |
_phi_jet[ multiplicity]->Fill( it->phi() );
|
| 396 |
_emfrac_jet[ multiplicity]->Fill( it->emEnergyFraction() );
|
| 397 |
_hadfrac_jet[multiplicity]->Fill( it->energyFractionHadronic() );
|
| 398 |
|
| 399 |
_n60_jet[ multiplicity]->Fill( it->n60() );
|
| 400 |
_n90_jet[ multiplicity]->Fill( it->n90() );
|
| 401 |
_area_jet[ multiplicity]->Fill( it->towersArea() );
|
| 402 |
|
| 403 |
FourierTransformation( multiplicity, *it, *EBRecHit );
|
| 404 |
}
|
| 405 |
++multiplicity;
|
| 406 |
ht += it->pt();
|
| 407 |
met += it->pt();//@@ just to kill warnings....
|
| 408 |
}
|
| 409 |
}
|
| 410 |
_jetmult->Fill( multiplicity );
|
| 411 |
_ht->Fill( ht );
|
| 412 |
if (multiplicity>1)
|
| 413 |
_dijet->Fill( sqrt( pow(patjets[0].energy()+patjets[1].energy(),2) -
|
| 414 |
pow(patjets[0].px()+patjets[1].px(),2) -
|
| 415 |
pow(patjets[0].py()+patjets[1].py(),2) -
|
| 416 |
pow(patjets[0].pz()+patjets[1].pz(),2)
|
| 417 |
) );
|
| 418 |
|
| 419 |
|
| 420 |
|
| 421 |
/*
|
| 422 |
//PAT MET
|
| 423 |
typedef vector<PatMet>::const_iterator PatJetCIter;
|
| 424 |
edm::Handle<vector<PatJet> > PatJets;
|
| 425 |
iEvent.getByLabel( _patJet, PatJets );
|
| 426 |
std::sort( PatJets->begin(), PatJets->end(), PtGreater());
|
| 427 |
PatJetCIter PatJetBegin = PatJets->begin();
|
| 428 |
PatJetCIter PatJetEnd = PatJets->end();
|
| 429 |
//MET
|
| 430 |
edm::Handle<reco::CaloMETCollection> recMet;
|
| 431 |
iEvent.getByLabel(_recMet,recMet);
|
| 432 |
|
| 433 |
for ( reco::CaloMETCollection::const_iterator it=recMet->begin();
|
| 434 |
it!=recMet->end(); ++it) {
|
| 435 |
met = it->pt();
|
| 436 |
_met->Fill( met );
|
| 437 |
_metx->Fill(met*sin(it->phi()));
|
| 438 |
_mety->Fill(met*cos(it->phi()));
|
| 439 |
break;
|
| 440 |
}
|
| 441 |
_metmult->Fill( recMet->size() );
|
| 442 |
_htmet->Fill( ht+met );
|
| 443 |
_metsig->Fill( met / (0.97*sqrt(ht)) );
|
| 444 |
if (multiplicity>1 ) {
|
| 445 |
_nv->Fill( (met*met-(CaloJets[0].pt()-CaloJets[1].pt())*(CaloJets[0].pt()-CaloJets[1].pt())) / (CaloJets[1].pt()*CaloJets[1].pt()) );
|
| 446 |
_nv2->Fill( (met*met-(CaloJets[0].pt()-CaloJets[1].pt())*(CaloJets[0].pt()-CaloJets[1].pt())) / (CaloJets[0].pt()*CaloJets[0].pt()) );
|
| 447 |
_dijet->Fill(sqrt(pow(CaloJets[0].energy()+CaloJets[1].energy(),2)
|
| 448 |
-(pow(CaloJets[0].px()+CaloJets[1].px(),2)
|
| 449 |
+pow(CaloJets[0].py()+CaloJets[1].py(),2)
|
| 450 |
+pow(CaloJets[0].pz()+CaloJets[1].pz(),2) ) ) );
|
| 451 |
}
|
| 452 |
*/
|
| 453 |
/*
|
| 454 |
//GenJets
|
| 455 |
multiplicity = 0;
|
| 456 |
met=0.0;
|
| 457 |
ht=0.0;
|
| 458 |
for (reco::GenJetCollection::const_iterator it=GenJets.begin();
|
| 459 |
it!=GenJets.end(); ++it) {
|
| 460 |
if (it->pt() > _jetminpt && fabs( it->eta() )<_jetmaxeta) {
|
| 461 |
if (multiplicity<_ngenjets ) {
|
| 462 |
_pt_genjet[ multiplicity]->Fill( it->pt() );
|
| 463 |
_eta_genjet[ multiplicity]->Fill( it->eta() );
|
| 464 |
_phi_genjet[ multiplicity]->Fill( it->phi() );
|
| 465 |
_emfrac_genjet[ multiplicity]->Fill( it->emEnergy() / it->energy() );
|
| 466 |
_hadfrac_genjet[ multiplicity]->Fill( it->hadEnergy() / it->energy() );
|
| 467 |
_invisible_genjet[multiplicity]->Fill( it->invisibleEnergy() / it->energy() );
|
| 468 |
_aux_genjet[ multiplicity]->Fill( it->auxiliaryEnergy() / it->energy() );
|
| 469 |
}
|
| 470 |
++multiplicity;
|
| 471 |
ht += it->pt();
|
| 472 |
}
|
| 473 |
}
|
| 474 |
_genjetmult->Fill( multiplicity );
|
| 475 |
|
| 476 |
//GenMET
|
| 477 |
edm::Handle<reco::GenMETCollection> genMet;
|
| 478 |
iEvent.getByLabel(_genMet,genMet);
|
| 479 |
|
| 480 |
for ( reco::GenMETCollection::const_iterator it=genMet->begin();
|
| 481 |
it!=genMet->end(); ++it) {
|
| 482 |
met = it->pt();
|
| 483 |
_genmet->Fill( met );
|
| 484 |
_genmetx->Fill(met*sin(it->phi()));
|
| 485 |
_genmety->Fill(met*cos(it->phi()));
|
| 486 |
break;
|
| 487 |
}
|
| 488 |
_genmetmult->Fill( recMet->size() );
|
| 489 |
_genht->Fill( ht );
|
| 490 |
_genhtmet->Fill( ht+met );
|
| 491 |
if (multiplicity>1 ) {
|
| 492 |
_gendijet->Fill(sqrt(pow(GenJets[0].energy()+GenJets[1].energy(),2)
|
| 493 |
-(pow(GenJets[0].px()+GenJets[1].px(),2)
|
| 494 |
+pow(GenJets[0].py()+GenJets[1].py(),2)
|
| 495 |
+pow(GenJets[0].pz()+GenJets[1].pz(),2) ) ) );
|
| 496 |
}
|
| 497 |
*/
|
| 498 |
|
| 499 |
|
| 500 |
|
| 501 |
}
|
| 502 |
|
| 503 |
|
| 504 |
// ------------ method called once each job just before starting event loop ------------
|
| 505 |
void
|
| 506 |
PATJetIDAnalyzer::beginJob(const edm::EventSetup&)
|
| 507 |
{
|
| 508 |
if( _hist.empty() )
|
| 509 |
return;
|
| 510 |
|
| 511 |
edm::Service<TFileService> fs;
|
| 512 |
if( !fs ){
|
| 513 |
throw edm::Exception( edm::errors::Configuration,
|
| 514 |
"TFile Service is not registered in cfg file" );
|
| 515 |
}
|
| 516 |
|
| 517 |
NameScheme hsusy("pat");
|
| 518 |
ofstream hist(_hist.c_str(), std::ios::out);
|
| 519 |
//CaloJets
|
| 520 |
for (unsigned i=0; i<_njets; ++i) {
|
| 521 |
_pt_jet[i] = fs->make<TH1F>(hsusy.name(hist, "ptj", i ), hsusy.name("PtJet", i), 100, 0. , 1000.);
|
| 522 |
_eta_jet[i] = fs->make<TH1F>(hsusy.name(hist, "etaj", i ), hsusy.name("EtaJet", i), 100, -5. , 5.);
|
| 523 |
_phi_jet[i] = fs->make<TH1F>(hsusy.name(hist, "phij", i ), hsusy.name("PhiJet", i), 70, -3.5, 3.5);
|
| 524 |
_emfrac_jet[i] = fs->make<TH1F>(hsusy.name(hist, "emfj", i ), hsusy.name("EmfJet", i), 50, 0. , 2.);
|
| 525 |
_hadfrac_jet[i] = fs->make<TH1F>(hsusy.name(hist, "hadfj",i ), hsusy.name("HadfJet",i), 50, 0. , 2.);
|
| 526 |
_n60_jet[i] = fs->make<TH1F>(hsusy.name(hist, "n60j", i ), hsusy.name("n60Jet", i), 50, 0. , 50.);
|
| 527 |
_n90_jet[i] = fs->make<TH1F>(hsusy.name(hist, "n90j", i ), hsusy.name("n90Jet", i), 50, 0. , 50.);
|
| 528 |
_area_jet[i] = fs->make<TH1F>(hsusy.name(hist, "areaj",i ), hsusy.name("AreaJet",i), 50, 0. , 5.);
|
| 529 |
|
| 530 |
_ft_energy[i] = fs->make<TH2F>(hsusy.name(hist, "ft_E_j",i ), hsusy.name("FtEJet",i), 16,-8,8,16,-8,8);
|
| 531 |
_ft_frequency[i]= fs->make<TH2F>(hsusy.name(hist, "ft_F_j",i ), hsusy.name("FtFJet",i), 16,-8,8,16,-8,8);
|
| 532 |
//_ft_frequency[i]= fs->make<TH1F>(hsusy.name(hist, "ft_F_j",i ), hsusy.name("FtFJet",i), 100,-50,50);
|
| 533 |
_ft_k[i] = fs->make<TH1F>(hsusy.name(hist, "ft_k_j",i ), hsusy.name("FtKJet",i), N_Fourier_Bins_1D,0.,0.5);
|
| 534 |
_ft_f[i] = fs->make<TH1F>(hsusy.name(hist, "ft_f_j",i ), hsusy.name("FtFJet",i), Nf_Fourier_Bins_1D+2,0.,N_Fourier_Bins_1D+1);
|
| 535 |
_ft_ksubavg[i] = fs->make<TH1F>(hsusy.name(hist, "ft_ksubavg_j",i ), hsusy.name("FtKsubavgJet",i), N_Fourier_Bins_1D,0.,0.5);
|
| 536 |
|
| 537 |
_noisecontrib[i]= fs->make<TH1F>(hsusy.name(hist, "noisecontrib_j",i ), hsusy.name("noisecontrib",i), 100,0.,20.);
|
| 538 |
_fto_n99[i] = fs->make<TH1I>(hsusy.name(hist, "fto_n99_j",i ), hsusy.name("fto_n99",i), 100,-1,500);
|
| 539 |
_fto_n95[i] = fs->make<TH1I>(hsusy.name(hist, "fto_n95_j",i ), hsusy.name("fto_n95",i), 100,-1,500);
|
| 540 |
_fto_n90[i] = fs->make<TH1I>(hsusy.name(hist, "fto_n90_j",i ), hsusy.name("fto_n90",i), 100,-1,500);
|
| 541 |
_fto_n60[i] = fs->make<TH1I>(hsusy.name(hist, "fto_n60_j",i ), hsusy.name("fto_n60",i), 100,-1,500);
|
| 542 |
_fto_n30[i] = fs->make<TH1I>(hsusy.name(hist, "fto_n30_j",i ), hsusy.name("fto_n30",i), 50,-1,500);
|
| 543 |
_fto_n10[i] = fs->make<TH1I>(hsusy.name(hist, "fto_n10_j",i ), hsusy.name("fto_n10",i), 50,-1,500);
|
| 544 |
_fto_F10[i] = fs->make<TH1F>(hsusy.name(hist, "fto_F10_j",i ), hsusy.name("fto_F10",i), 100,0.,2000.);
|
| 545 |
_fto_Fs10[i] = fs->make<TH1F>(hsusy.name(hist, "fto_Fs10_j",i),hsusy.name("fto_Fs10",i), 100,0.,100.);
|
| 546 |
_fto_LowFvsHiF[i]=fs->make<TH1F>(hsusy.name(hist, "fto_LowFvsHiF_j",i ), hsusy.name("LowFvsHiF",i), 100,-5.,5.);
|
| 547 |
_fto_LowFovHiF[i]=fs->make<TH1F>(hsusy.name(hist, "fto_LowFovHiF_j",i ), hsusy.name("LowFovHiF",i), 100,-5.,5.);
|
| 548 |
_fto_det16[i] =fs->make<TH1F>(hsusy.name(hist, "fto_Det16_j",i ), hsusy.name("Det16",i), 100,-1.0E08,1.0E08);
|
| 549 |
_fto_det8[i] =fs->make<TH1F>(hsusy.name(hist, "fto_Det8_j",i ), hsusy.name("Det8",i), 100,-1.0E08,1.0E08);
|
| 550 |
|
| 551 |
}
|
| 552 |
_jetmult = fs->make<TH1F>(hsusy.name(hist, "multj" ), hsusy.name("JetMultiplicity"), 21, -0.5, 20.5);
|
| 553 |
_ht = fs->make<TH1F>(hsusy.name(hist, "ht" ), hsusy.name("Ht"), 100, 0.0, 2000.);
|
| 554 |
_htmet = fs->make<TH1F>(hsusy.name(hist, "htmet" ), hsusy.name("HtMet"), 100, 0.0, 2000.);
|
| 555 |
_dijet = fs->make<TH1F>(hsusy.name(hist, "dijet" ), hsusy.name("DiJet"), 100, 0.0, 2000.);
|
| 556 |
//GenJets
|
| 557 |
for (unsigned i=0; i<_ngenjets; ++i) {
|
| 558 |
_pt_genjet[i] = fs->make<TH1F>(hsusy.name(hist, "ptgj", i ), hsusy.name("PtGenJet", i), 100, 0. , 1000.);
|
| 559 |
_eta_genjet[i] = fs->make<TH1F>(hsusy.name(hist, "etagj", i ), hsusy.name("EtaGenJet", i), 100, -5. , 5.);
|
| 560 |
_phi_genjet[i] = fs->make<TH1F>(hsusy.name(hist, "phigj", i ), hsusy.name("PhiGenJet", i), 70, -3.5, 3.5);
|
| 561 |
_emfrac_genjet[i] = fs->make<TH1F>(hsusy.name(hist, "emfgj", i ), hsusy.name("EmfGenJet", i), 50, 0. , 5.);
|
| 562 |
_hadfrac_genjet[i] = fs->make<TH1F>(hsusy.name(hist, "hadfgj",i ), hsusy.name("HadfGenJet",i), 50, 0. , 5.);
|
| 563 |
_invisible_genjet[i]=fs->make<TH1F>(hsusy.name(hist, "invisgj",i), hsusy.name("InvisiblefGenJet", i), 50, 0.,5.);
|
| 564 |
_aux_genjet[i] = fs->make<TH1F>(hsusy.name(hist, "auxgj", i ), hsusy.name("AuxfGenJet", i), 50, 0. , 5.);
|
| 565 |
}
|
| 566 |
_genjetmult = fs->make<TH1F>(hsusy.name(hist, "multgj" ), hsusy.name("GenJetMultiplicity"), 21, -0.5, 20.5);
|
| 567 |
_genht = fs->make<TH1F>(hsusy.name(hist, "htgj" ), hsusy.name("GenHt"), 100, 0.0, 2000.);
|
| 568 |
_genhtmet = fs->make<TH1F>(hsusy.name(hist, "htmetgj"), hsusy.name("GenHtMet"), 100, 0.0, 2000.);
|
| 569 |
_gendijet = fs->make<TH1F>(hsusy.name(hist, "dijetgj"), hsusy.name("GenDiJet"), 100, 0.0, 2000.);
|
| 570 |
//MET
|
| 571 |
_met = fs->make<TH1F>(hsusy.name(hist, "met" ), hsusy.name("Met"), 100, 0.0, 1000.);
|
| 572 |
_metmult = fs->make<TH1F>(hsusy.name(hist, "metmult"),hsusy.name("MetMult"), 11, -0.5, 10.5);
|
| 573 |
_metx = fs->make<TH1F>(hsusy.name(hist, "metx" ), hsusy.name("MetX"), 100, -1000.0, 1000.);
|
| 574 |
_mety = fs->make<TH1F>(hsusy.name(hist, "mety" ), hsusy.name("MetY"), 100, -1000.0, 1000.);
|
| 575 |
_metsig = fs->make<TH1F>(hsusy.name(hist, "metsig"), hsusy.name("MetSig"), 100, 0.0, 100.);
|
| 576 |
//GenMET
|
| 577 |
_genmet = fs->make<TH1F>(hsusy.name(hist, "metgj" ), hsusy.name("GenMet"), 100, 0.0, 1000.);
|
| 578 |
_genmetmult = fs->make<TH1F>(hsusy.name(hist, "metmultgj"),hsusy.name("GenMetMult"), 11, -0.5, 10.5);
|
| 579 |
_genmetx = fs->make<TH1F>(hsusy.name(hist, "metxgj" ), hsusy.name("GenMetX"), 100, -1000.0, 1000.);
|
| 580 |
_genmety = fs->make<TH1F>(hsusy.name(hist, "metygj" ), hsusy.name("GenMetY"), 100, -1000.0, 1000.);
|
| 581 |
//other
|
| 582 |
_nv = fs->make<TH1F>(hsusy.name(hist, "nv" ), hsusy.name("NV"), 100, -5.0, 5.);
|
| 583 |
_nv2 = fs->make<TH1F>(hsusy.name(hist, "nv2" ), hsusy.name("NV2"), 100, -5.0, 5.);
|
| 584 |
//matching
|
| 585 |
_GenOnCalo_match = fs->make<TH2F>(hsusy.name(hist,"gencalo"), hsusy.name("GenOnCalo"), 60, 0., 3., 120, -1, 5.);
|
| 586 |
_CaloOnGen_match = fs->make<TH2F>(hsusy.name(hist,"calogen"), hsusy.name("CaloOnGen"), 60, 0., 3., 120, -1, 5.);
|
| 587 |
_GenVsMatched_pt = fs->make<TH2F>(hsusy.name(hist,"genmatched_pt"), hsusy.name("GenVsMatched_Pt") , 100, 0. , 2000. , 100, 0. , 2000.);
|
| 588 |
_GenVsMatched_eta = fs->make<TH2F>(hsusy.name(hist,"genmatched_eta"), hsusy.name("GenVsMatched_Eta"), 100, -5. , 5. , 100, -5. , 5.);
|
| 589 |
_GenVsMatched_phi = fs->make<TH2F>(hsusy.name(hist,"genmatched_phi"), hsusy.name("GenVsMatched_Phi"), 70, -3.5, 3.5, 70, -3.5, 3.5);
|
| 590 |
_RecoEff_pt = fs->make<TH1F>(hsusy.name(hist, "recoeffpt") , hsusy.name("RecoEff_Pt") , 100, 0. , 2000. );
|
| 591 |
_RecoEff_eta = fs->make<TH1F>(hsusy.name(hist, "recoeffeta"), hsusy.name("RecoEff_Eta"), 100, -5. , 5. );
|
| 592 |
_RecoEff_phi = fs->make<TH1F>(hsusy.name(hist, "recoeffphi"), hsusy.name("RecoEff_Phi"), 70, -3.5, 3.5);
|
| 593 |
//helper histograms
|
| 594 |
_RecoTot_pt = new TH1F("recototpt" , "RecoTot_Pt" , 100, 0. , 2000. );
|
| 595 |
_RecoTot_eta= new TH1F("recototeta", "RecoTot_Eta", 100, -5. , 5. );
|
| 596 |
_RecoTot_phi= new TH1F("recototphi", "RecoTot_Phi", 70, -3.5, 3.5);
|
| 597 |
}
|
| 598 |
|
| 599 |
// ------------ method called once each job just after ending the event loop ------------
|
| 600 |
void
|
| 601 |
PATJetIDAnalyzer::endJob() {
|
| 602 |
_RecoTot_pt->Sumw2();
|
| 603 |
_RecoEff_pt->Sumw2();
|
| 604 |
_RecoEff_pt->Divide(_RecoTot_pt);
|
| 605 |
_RecoTot_eta->Sumw2();
|
| 606 |
_RecoEff_eta->Sumw2();
|
| 607 |
_RecoEff_eta->Divide(_RecoTot_eta);
|
| 608 |
_RecoTot_phi->Sumw2();
|
| 609 |
_RecoEff_phi->Sumw2();
|
| 610 |
_RecoEff_phi->Divide(_RecoTot_phi);
|
| 611 |
}
|
| 612 |
|
| 613 |
// ------------ redo matching
|
| 614 |
void
|
| 615 |
PATJetIDAnalyzer::makeMatchingMaps(edm::Handle<reco::GenJetCollection> GenJets, edm::Handle<reco::CaloJetCollection> CaloJets)
|
| 616 |
{
|
| 617 |
MatchingMapGen.clear();
|
| 618 |
MatchingMapJet.clear();
|
| 619 |
|
| 620 |
// Loop over generator Jets
|
| 621 |
for(reco::GenJetCollection::const_iterator genjet = GenJets->begin(); genjet!=GenJets->end(); ++genjet) {
|
| 622 |
const reco::CaloJet* jetbest=0;
|
| 623 |
const reco::CaloJet* jetmatched=0;
|
| 624 |
double minDelR= 0.3;
|
| 625 |
double minDelE=-1/4.;
|
| 626 |
double maxDelE= 3.0;
|
| 627 |
double mindE=999.;
|
| 628 |
double mindR=999.;
|
| 629 |
double delR;
|
| 630 |
double delE;
|
| 631 |
double E1=(*genjet).energy();
|
| 632 |
bool matched = false;
|
| 633 |
// Loop over Calo jets
|
| 634 |
for(reco::CaloJetCollection::const_iterator calojet = CaloJets->begin(); calojet!=CaloJets->end(); ++calojet){
|
| 635 |
double E2=(*calojet).energy();
|
| 636 |
delR=deltaR(*calojet,*genjet);
|
| 637 |
delE=(E2-E1)/E1;
|
| 638 |
// This is the matching condition
|
| 639 |
if(delR<minDelR && delE>minDelE && delE<maxDelE){
|
| 640 |
matched=true;
|
| 641 |
minDelR=delR;
|
| 642 |
jetmatched=&(*calojet);
|
| 643 |
}
|
| 644 |
if(delR<mindR){
|
| 645 |
mindR=delR;
|
| 646 |
mindE=delE;
|
| 647 |
jetbest=&(*calojet);
|
| 648 |
}
|
| 649 |
}
|
| 650 |
if (matched) {
|
| 651 |
MatchingMapGen.insert(make_pair(&(*genjet),&(*jetmatched)));
|
| 652 |
_GenVsMatched_pt->Fill((*genjet).pt(),(*jetmatched).pt());
|
| 653 |
_GenVsMatched_eta->Fill((*genjet).eta(),(*jetmatched).eta());
|
| 654 |
_GenVsMatched_phi->Fill((*genjet).phi(),(*jetmatched).phi());
|
| 655 |
_RecoEff_pt->Fill((*genjet).pt());
|
| 656 |
_RecoEff_eta->Fill((*genjet).eta());
|
| 657 |
_RecoEff_phi->Fill((*genjet).phi());
|
| 658 |
}
|
| 659 |
_RecoTot_pt->Fill((*genjet).pt());
|
| 660 |
_RecoTot_eta->Fill((*genjet).eta());
|
| 661 |
_RecoTot_phi->Fill((*genjet).phi());
|
| 662 |
_GenOnCalo_match->Fill(mindR,mindE);
|
| 663 |
}
|
| 664 |
|
| 665 |
// Loop over Calo Jets
|
| 666 |
for(reco::CaloJetCollection::const_iterator calojet = CaloJets->begin(); calojet!=CaloJets->end(); ++calojet) {
|
| 667 |
const reco::GenJet* jetbest=0;
|
| 668 |
const reco::GenJet* jetmatched=0;
|
| 669 |
double minDelR= 0.3;
|
| 670 |
double minDelE=-1/4.;
|
| 671 |
double maxDelE= 3.0;
|
| 672 |
double mindE=999.;
|
| 673 |
double mindR=999.;
|
| 674 |
double delR;
|
| 675 |
double delE;
|
| 676 |
double E1=(*calojet).energy();
|
| 677 |
bool matched = false;
|
| 678 |
// Loop over jets
|
| 679 |
for(reco::GenJetCollection::const_iterator genjet = GenJets->begin(); genjet!=GenJets->end(); ++genjet){
|
| 680 |
double E2=(*genjet).energy();
|
| 681 |
delR=deltaR(*genjet,*calojet);
|
| 682 |
delE=(E2-E1)/E1;
|
| 683 |
// This is the matching condition
|
| 684 |
if(delR<minDelR && delE>minDelE && delE<maxDelE){
|
| 685 |
matched=true;
|
| 686 |
minDelR=delR;
|
| 687 |
jetmatched=&(*genjet);
|
| 688 |
}
|
| 689 |
if(delR<mindR){
|
| 690 |
mindR=delR;
|
| 691 |
mindE=delE;
|
| 692 |
jetbest=&(*genjet);
|
| 693 |
}
|
| 694 |
}
|
| 695 |
if (matched) {
|
| 696 |
MatchingMapJet.insert(make_pair(&(*calojet),&(*jetmatched)));
|
| 697 |
}
|
| 698 |
_CaloOnGen_match->Fill(mindR,mindE);
|
| 699 |
}
|
| 700 |
|
| 701 |
}
|
| 702 |
|
| 703 |
//define this as a plug-in
|
| 704 |
DEFINE_FWK_MODULE(PATJetIDAnalyzer);
|
