30 jetsToken( consumes<edm::
View<
reco::
Jet>>( iConfig.getParameter<edm::
InputTag>(
"srcJets"))),
33 rhoToken( consumes<double>( iConfig.getParameter<edm::
InputTag>(
"srcRho"))),
35 systLabel( iConfig.getParameter<std::
string>(
"systematicsLabel")),
36 useQC( iConfig.getParameter<bool>(
"useQualityCuts")),
38 produceSyst( systLabel !=
"")
40 produces<edm::ValueMap<float>>(
"qgLikelihood");
41 produces<edm::ValueMap<float>>(
"axis2");
42 produces<edm::ValueMap<int>>(
"mult");
43 produces<edm::ValueMap<float>>(
"ptD");
45 produces<edm::ValueMap<float>>(
"qgLikelihoodSmearedQuark");
46 produces<edm::ValueMap<float>>(
"qgLikelihoodSmearedGluon");
47 produces<edm::ValueMap<float>>(
"qgLikelihoodSmearedAll");
56 std::vector<float>* qgProduct =
new std::vector<float>;
57 std::vector<float>* axis2Product =
new std::vector<float>;
58 std::vector<int>* multProduct =
new std::vector<int>;
59 std::vector<float>* ptDProduct =
new std::vector<float>;
60 std::vector<float>* smearedQuarkProduct =
new std::vector<float>;
61 std::vector<float>* smearedGluonProduct =
new std::vector<float>;
62 std::vector<float>* smearedAllProduct =
new std::vector<float>;
79 for(
auto jet = jets->begin();
jet != jets->end(); ++
jet){
82 float ptD, axis2;
int mult;
89 qgProduct->push_back(qgValue);
95 axis2Product->push_back(axis2);
96 multProduct->push_back(mult);
97 ptDProduct->push_back(ptD);
100 putInEvent(
"qgLikelihood", jets, qgProduct, iEvent);
101 putInEvent(
"axis2", jets, axis2Product, iEvent);
102 putInEvent(
"mult", jets, multProduct, iEvent);
105 putInEvent(
"qgLikelihoodSmearedQuark", jets, smearedQuarkProduct, iEvent);
106 putInEvent(
"qgLikelihoodSmearedGluon", jets, smearedGluonProduct, iEvent);
107 putInEvent(
"qgLikelihoodSmearedAll", jets, smearedAllProduct, iEvent);
115 filler.
insert(
jets, product->begin(), product->end());
117 iEvent.
put(out, name);
127 else throw cms::Exception(
"WrongJetCollection",
"Jet constituents are not particle flow candidates");
136 float sum_weight = 0., sum_deta = 0., sum_dphi = 0., sum_deta2 = 0., sum_dphi2 = 0., sum_detadphi = 0., sum_pt = 0.;
145 if(!(
part->fromPV() > 1 &&
part->trackHighPurity()))
continue;
147 if((
part->dz()*
part->dz())/(
part->dzError()*
part->dzError()) > 25.)
continue;
151 if(
part->pt() < 1.0)
continue;
159 auto vtxLead = vC->begin();
160 auto vtxClose = vC->begin();
161 for(
auto vtx = vC->begin(); vtx != vC->end(); ++vtx){
162 if(fabs(itrk->dz(vtx->position())) < fabs(itrk->dz(vtxClose->position()))) vtxClose = vtx;
167 float dz = itrk->dz(vtxClose->position());
168 float d0 = itrk->dxy(vtxClose->position());
169 float dz_sigma_square =
pow(itrk->dzError(),2) +
pow(vtxClose->zError(),2);
170 float d0_sigma_square =
pow(itrk->d0Error(),2) +
pow(vtxClose->xError(),2) +
pow(vtxClose->yError(),2);
171 if(dz*dz/dz_sigma_square > 25.)
continue;
172 if(d0*d0/d0_sigma_square < 25.) ++
mult;
175 if(
part->pt() < 1.0)
continue;
180 float deta = daughter->eta() - jet->
eta();
182 float partPt = daughter->pt();
183 float weight = partPt*partPt;
189 sum_deta2 += deta*deta*
weight;
190 sum_detadphi += deta*dphi*
weight;
191 sum_dphi2 += dphi*dphi*
weight;
195 float a = 0.,
b = 0.,
c = 0.;
196 float ave_deta = 0., ave_dphi = 0., ave_deta2 = 0., ave_dphi2 = 0.;
198 ave_deta = sum_deta/sum_weight;
199 ave_dphi = sum_dphi/sum_weight;
200 ave_deta2 = sum_deta2/sum_weight;
201 ave_dphi2 = sum_dphi2/sum_weight;
202 a = ave_deta2 - ave_deta*ave_deta;
203 b = ave_dphi2 - ave_dphi*ave_dphi;
204 c = -(sum_detadphi/sum_weight - ave_deta*ave_dphi);
207 float axis2 = (a+
b-delta > 0 ?
sqrt(0.5*(a+
b-delta)) : 0);
208 float ptD = (sum_weight > 0 ?
sqrt(sum_weight)/sum_pt : 0);
209 return std::make_tuple(mult, ptD, axis2);
220 desc.
add<
bool>(
"useQualityCuts");
222 desc.
add<
edm::InputTag>(
"srcVertexCollection")->setComment(
"Ignored for miniAOD, possible to keep empty");
223 descriptions.
add(
"QGTagger", desc);
bool isNonnull() const
Checks for non-null.
bool getByToken(EDGetToken token, Handle< PROD > &result) const
#define DEFINE_FWK_MODULE(type)
Base class for all types of Jets.
void insert(const H &h, I begin, I end)
edm::EDGetTokenT< reco::VertexCollection > vertexToken
std::vector< Vertex > VertexCollection
collection of Vertex objects
virtual void produce(edm::Event &, const edm::EventSetup &)
Produce qgLikelihood using {mult, ptD, -log(axis2)}.
QGLikelihoodCalculator * qgLikelihood
QGTagger(const edm::ParameterSet &)
virtual double eta() const
momentum pseudorapidity
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
float computeQGLikelihood(edm::ESHandle< QGLikelihoodObject > &QGLParamsColl, float pt, float eta, float rho, std::vector< float > vars)
Compute likelihood for a jet using the QGLikelihoodObject information and a set of variables...
void get(HolderT &iHolder) const
bool isPackedCandidate(const reco::Candidate *candidate)
Function to tell us if we are using packedCandidates, only test for first candidate.
ParameterDescriptionBase * add(U const &iLabel, T const &value)
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Descriptions method.
std::tuple< int, float, float > calcVariables(const reco::Jet *, edm::Handle< reco::VertexCollection > &)
Calculation of axis2, mult and ptD.
edm::EDGetTokenT< edm::View< reco::Jet > > jetsToken
double deltaPhi(double phi1, double phi2)
static TrackQuality qualityByName(const std::string &name)
void add(std::string const &label, ParameterSetDescription const &psetDescription)
edm::EDGetTokenT< double > rhoToken
bool weAreUsingPackedCandidates
Particle reconstructed by the particle flow algorithm.
void putInEvent(std::string, const edm::Handle< edm::View< reco::Jet >> &, std::vector< T > *, edm::Event &)
Function to put product into event.
bool weStillNeedToCheckJetCandidates
edm::EDGetTokenT< reco::JetCorrector > jetCorrectorToken
virtual std::vector< const reco::Candidate * > getJetConstituentsQuick() const
quick list of constituents
virtual double phi() const
momentum azimuthal angle
Power< A, B >::type pow(const A &a, const B &b)
float systematicSmearing(edm::ESHandle< QGLikelihoodSystematicsObject > &QGLParamsColl, float pt, float eta, float rho, float qgValue, int qgIndex)