239 for(
unsigned i=0;
i<
pfJets.size();
i++) {
242 unsigned highJets = 0;
243 for(
unsigned j=0; j<
pfJets.size(); j++) {
250 double rec_pt = pfj.
pt();
251 double rec_eta = pfj.
eta();
252 double rec_phi = pfj.
phi();
269 bool Forward =
false;
270 if (
std::abs(rec_eta) < 1.4 ) Barrel =
true;
288 double true_E = truth->
p();
289 double true_pt = truth->
pt();
290 double true_eta = truth->
eta();
291 double true_phi = truth->
phi();
294 else {pt_denom = true_pt;}
296 double true_ChargedHadEnergy;
297 double true_NeutralHadEnergy;
298 double true_NeutralEmEnergy;
299 gettrue (truth, true_ChargedHadEnergy, true_NeutralHadEnergy, true_NeutralEmEnergy);
300 double true_NeutralEnergy = true_NeutralHadEnergy + true_NeutralEmEnergy;
304 double rec_NeutralEnergy = rec_NeutralHadEnergy + rec_NeutralEmEnergy;
307 std::vector <unsigned int> chMult(9, static_cast<unsigned int>(0));
308 for (
unsigned ic = 0; ic < constituents.size (); ++ic) {
309 if ( constituents[ic]->particleId() > 3 )
continue;
311 if ( trackRef.
isNull() ) {
317 unsigned int iter = trackRef->algo()>6 ? 6: trackRef->algo();
328 double cut1 = 0.0001;
329 double cut2 = 0.0001;
330 double cut3 = 0.0001;
331 double cut4 = 0.0001;
332 double cut5 = 0.0001;
333 double cut6 = 0.0001;
334 double cut7 = 0.0001;
336 double resChargedHadEnergy= 0.;
337 double resNeutralHadEnergy= 0.;
338 double resNeutralEmEnergy= 0.;
339 double resNeutralEnergy= 0.;
341 double resHCALEnergy = 0.;
342 double resNEUTEnergy = 0.;
343 if ( rec_NeutralHadEnergy > cut6 && rec_ChargedHadEnergy < cut1 ) {
344 double true_NEUTEnergy = true_NeutralHadEnergy + true_NeutralEmEnergy;
345 double true_HCALEnergy = true_NEUTEnergy - rec_NeutralEmEnergy;
346 double rec_NEUTEnergy = rec_NeutralHadEnergy+rec_NeutralEmEnergy;
347 double rec_HCALEnergy = rec_NeutralHadEnergy;
348 resHCALEnergy = (rec_HCALEnergy-true_HCALEnergy)/rec_HCALEnergy;
349 resNEUTEnergy = (rec_NEUTEnergy-true_NEUTEnergy)/rec_NEUTEnergy;
350 if ( rec_NeutralEmEnergy > cut7 ) {
358 if (true_pt > 0.0001){
359 resPt = (rec_pt -true_pt)/true_pt ;
361 if (true_ChargedHadEnergy > cut1){
362 resChargedHadEnergy = (rec_ChargedHadEnergy- true_ChargedHadEnergy)/true_ChargedHadEnergy;
364 if (true_NeutralHadEnergy > cut2){
365 resNeutralHadEnergy = (rec_NeutralHadEnergy- true_NeutralHadEnergy)/true_NeutralHadEnergy;
368 if (rec_NeutralHadEnergy > cut3){
369 resNeutralHadEnergy = (rec_NeutralHadEnergy- true_NeutralHadEnergy)/rec_NeutralHadEnergy;
371 if (true_NeutralEmEnergy > cut4){
372 resNeutralEmEnergy = (rec_NeutralEmEnergy- true_NeutralEmEnergy)/true_NeutralEmEnergy;
374 if (true_NeutralEnergy > cut5){
375 resNeutralEnergy = (rec_NeutralEnergy- true_NeutralEnergy)/true_NeutralEnergy;
382 if ( ( resPt > 0.2 && true_pt > 100. ) ||
383 ( resPt < -0.5 && true_pt > 100. ) ) {
388 <<
" resPt = " << resPt
389 <<
" resCharged " << resChargedHadEnergy
390 <<
" resNeutralHad " << resNeutralHadEnergy
391 <<
" resNeutralEm " << resNeutralEmEnergy
392 <<
" pT (T/R) " << true_pt <<
"/" << rec_pt
393 <<
" Eta (T/R) " << truth->
eta() <<
"/" << rec_eta
394 <<
" Phi (T/R) " << truth->
phi() <<
"/" << rec_phi
400 for(
unsigned j=0; j<
pfJets.size(); j++) {
402 if ( j !=
i &&
algo_->
deltaR(&pfj,&pfo) < dRo && pfo.
pt() > 0.25*pfj.
pt()) {
412 for(
unsigned j=0; j<
genJets.size(); j++) {
414 if ( gjo != truth &&
algo_->
deltaR(truth,gjo) < dRgo && gjo->
pt() > 0.25*truth->
pt() ) {
421 std::cout <<
"Excess probably due to overlapping jets (DR = " <<
algo_->
deltaR(genoj,pfoj) <<
")," 422 <<
" at DeltaR(T/R) = " << dRgo <<
"/" << dRo
423 <<
" with pT(T/R) " << genoj->
pt() <<
"/" << pfoj->
pt()
424 <<
" and Eta (T/R) " << genoj->
eta() <<
"/" << pfoj->
eta()
425 <<
" and Phi (T/R) " << genoj->
phi() <<
"/" << pfoj->
phi()
434 cout <<
i <<
" =========PFJet Pt "<< rec_pt
435 <<
" eta " << rec_eta
436 <<
" phi " << rec_phi
437 <<
" Charged Had Energy " << rec_ChargedHadEnergy
438 <<
" Neutral Had Energy " << rec_NeutralHadEnergy
439 <<
" Neutral elm Energy " << rec_NeutralEmEnergy << endl;
440 cout <<
" matching Gen Jet Pt " << true_pt
441 <<
" eta " << truth->
eta()
442 <<
" phi " << truth->
phi()
443 <<
" Charged Had Energy " << true_ChargedHadEnergy
444 <<
" Neutral Had Energy " << true_NeutralHadEnergy
445 <<
" Neutral elm Energy " << true_NeutralEmEnergy << endl;
451 cout <<
"==============deltaR " << deltaR <<
" resPt " << resPt
452 <<
" resChargedHadEnergy " << resChargedHadEnergy
453 <<
" resNeutralHadEnergy " << resNeutralHadEnergy
454 <<
" resNeutralEmEnergy " << resNeutralEmEnergy
467 hNCHvsEta->Fill(true_eta, rec_ChargedMultiplicity);
477 if(plot2)
hRCHEvsEta->Fill(true_eta, resChargedHadEnergy);
478 if(plot5)
hRNeutvsEta->Fill(true_eta, resNeutralEnergy);
493 if ( pt_denom > 20. && pt_denom < 40. )
hBRPt20_40->Fill (resPt);
494 if ( pt_denom > 40. && pt_denom < 60. )
hBRPt40_60->Fill (resPt);
495 if ( pt_denom > 60. && pt_denom < 80. )
hBRPt60_80->Fill (resPt);
496 if ( pt_denom > 80. && pt_denom < 100. )
hBRPt80_100->Fill (resPt);
497 if ( pt_denom > 100. && pt_denom < 150. )
hBRPt100_150->Fill (resPt);
498 if ( pt_denom > 150. && pt_denom < 200. )
hBRPt150_200->Fill (resPt);
499 if ( pt_denom > 200. && pt_denom < 250. )
hBRPt200_250->Fill (resPt);
500 if ( pt_denom > 250. && pt_denom < 300. )
hBRPt250_300->Fill (resPt);
501 if ( pt_denom > 300. && pt_denom < 400. )
hBRPt300_400->Fill (resPt);
502 if ( pt_denom > 400. && pt_denom < 500. )
hBRPt400_500->Fill (resPt);
503 if ( pt_denom > 500. && pt_denom < 750. )
hBRPt500_750->Fill (resPt);
504 if ( pt_denom > 750. && pt_denom < 1250. )
hBRPt750_1250->Fill (resPt);
505 if ( pt_denom > 1250. && pt_denom < 2000. )
hBRPt1250_2000->Fill (resPt);
506 if ( pt_denom > 2000. && pt_denom < 5000. )
hBRPt2000_5000->Fill (resPt);
507 hBNCH->Fill(rec_ChargedMultiplicity);
516 hBNCHvsPt->Fill(pt_denom,rec_ChargedMultiplicity);
523 if(plot2)
hBRCHE->Fill(resChargedHadEnergy);
524 if(plot3)
hBRNHE->Fill(resNeutralHadEnergy);
525 if(plot4)
hBRNEE->Fill(resNeutralEmEnergy);
526 if(plot5)
hBRneut->Fill(resNeutralEnergy);
527 if(plot1)
hBRPtvsPt->Fill(pt_denom, resPt);
528 if(plot2)
hBRCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
529 if(plot3)
hBRNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
530 if(plot4)
hBRNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
531 if(plot5)
hBRneutvsPt->Fill(pt_denom, resNeutralEnergy);
546 if ( pt_denom > 20. && pt_denom < 40. )
hERPt20_40->Fill (resPt);
547 if ( pt_denom > 40. && pt_denom < 60. )
hERPt40_60->Fill (resPt);
548 if ( pt_denom > 60. && pt_denom < 80. )
hERPt60_80->Fill (resPt);
549 if ( pt_denom > 80. && pt_denom < 100. )
hERPt80_100->Fill (resPt);
550 if ( pt_denom > 100. && pt_denom < 150. )
hERPt100_150->Fill (resPt);
551 if ( pt_denom > 150. && pt_denom < 200. )
hERPt150_200->Fill (resPt);
552 if ( pt_denom > 200. && pt_denom < 250. )
hERPt200_250->Fill (resPt);
553 if ( pt_denom > 250. && pt_denom < 300. )
hERPt250_300->Fill (resPt);
554 if ( pt_denom > 300. && pt_denom < 400. )
hERPt300_400->Fill (resPt);
555 if ( pt_denom > 400. && pt_denom < 500. )
hERPt400_500->Fill (resPt);
556 if ( pt_denom > 500. && pt_denom < 750. )
hERPt500_750->Fill (resPt);
557 if ( pt_denom > 750. && pt_denom < 1250. )
hERPt750_1250->Fill (resPt);
558 if ( pt_denom > 1250. && pt_denom < 2000. )
hERPt1250_2000->Fill (resPt);
559 if ( pt_denom > 2000. && pt_denom < 5000. )
hERPt2000_5000->Fill (resPt);
560 hENCH->Fill(rec_ChargedMultiplicity);
561 hENCHvsPt->Fill(pt_denom,rec_ChargedMultiplicity);
576 if(plot2)
hERCHE->Fill(resChargedHadEnergy);
577 if(plot3)
hERNHE->Fill(resNeutralHadEnergy);
578 if(plot4)
hERNEE->Fill(resNeutralEmEnergy);
579 if(plot5)
hERneut->Fill(resNeutralEnergy);
580 if(plot1)
hERPtvsPt->Fill(pt_denom, resPt);
581 if(plot2)
hERCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
582 if(plot3)
hERNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
583 if(plot4)
hERNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
584 if(plot5)
hERneutvsPt->Fill(pt_denom, resNeutralEnergy);
598 if ( pt_denom > 20. && pt_denom < 40. )
hFRPt20_40->Fill (resPt);
599 if ( pt_denom > 40. && pt_denom < 60. )
hFRPt40_60->Fill (resPt);
600 if ( pt_denom > 60. && pt_denom < 80. )
hFRPt60_80->Fill (resPt);
601 if ( pt_denom > 80. && pt_denom < 100. )
hFRPt80_100->Fill (resPt);
602 if ( pt_denom > 100. && pt_denom < 150. )
hFRPt100_150->Fill (resPt);
603 if ( pt_denom > 150. && pt_denom < 200. )
hFRPt150_200->Fill (resPt);
604 if ( pt_denom > 200. && pt_denom < 250. )
hFRPt200_250->Fill (resPt);
605 if ( pt_denom > 250. && pt_denom < 300. )
hFRPt250_300->Fill (resPt);
606 if ( pt_denom > 300. && pt_denom < 400. )
hFRPt300_400->Fill (resPt);
607 if ( pt_denom > 400. && pt_denom < 500. )
hFRPt400_500->Fill (resPt);
608 if ( pt_denom > 500. && pt_denom < 750. )
hFRPt500_750->Fill (resPt);
609 if ( pt_denom > 750. && pt_denom < 1250. )
hFRPt750_1250->Fill (resPt);
610 if ( pt_denom > 1250. && pt_denom < 2000. )
hFRPt1250_2000->Fill (resPt);
611 if ( pt_denom > 2000. && pt_denom < 5000. )
hFRPt2000_5000->Fill (resPt);
618 if(plot2)
hFRCHE->Fill(resChargedHadEnergy);
619 if(plot3)
hFRNHE->Fill(resNeutralHadEnergy);
620 if(plot4)
hFRNEE->Fill(resNeutralEmEnergy);
621 if(plot5)
hFRneut->Fill(resNeutralEnergy);
622 if(plot1)
hFRPtvsPt->Fill(pt_denom, resPt);
623 if(plot2)
hFRCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
624 if(plot3)
hFRNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
625 if(plot4)
hFRNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
626 if(plot5)
hFRneutvsPt->Fill(pt_denom, resNeutralEnergy);
double eta() const final
momentum pseudorapidity
double pt() const final
transverse momentum
int chargedMultiplicity() const
chargedMultiplicity
Jets made from PFObjects.
float neutralEmEnergy() const
neutralEmEnergy
void gettrue(const reco::GenJet *truth, double &true_ChargedHadEnergy, double &true_NeutralHadEnergy, double &true_NeutralEmEnergy)
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
static double deltaR(const T *, const U *)
Abs< T >::type abs(const T &t)
Jets made from MC generator particles.
bool isNull() const
Checks for null.
double deltaR(double eta1, double eta2, double phi1, double phi2)
double p() const final
magnitude of momentum vector
double resChargedHadEnergyMax_
static const Collection::value_type * matchByDeltaR(const T *, const Collection *)
double resNeutralHadEnergyMax_
void printPFJet(const reco::PFJet *)
virtual std::vector< reco::PFCandidatePtr > getPFConstituents() const
get all constituents
float neutralHadronEnergy() const
neutralHadronEnergy
double resNeutralEmEnergyMax_
double phi() const final
momentum azimuthal angle
void printGenJet(const reco::GenJet *)
float chargedHadronEnergy() const
chargedHadronEnergy