238 for(
unsigned i=0;
i<
pfJets.size();
i++) {
241 unsigned highJets = 0;
242 for(
unsigned j=0;
j<
pfJets.size();
j++) {
249 double rec_pt = pfj.
pt();
250 double rec_eta = pfj.
eta();
251 double rec_phi = pfj.
phi();
268 bool Forward =
false;
269 if (
std::abs(rec_eta) < 1.4 ) Barrel =
true;
287 double true_E = truth->
p();
288 double true_pt = truth->
pt();
289 double true_eta = truth->
eta();
290 double true_phi = truth->
phi();
293 else {pt_denom = true_pt;}
295 double true_ChargedHadEnergy;
296 double true_NeutralHadEnergy;
297 double true_NeutralEmEnergy;
298 gettrue (truth, true_ChargedHadEnergy, true_NeutralHadEnergy, true_NeutralEmEnergy);
299 double true_NeutralEnergy = true_NeutralHadEnergy + true_NeutralEmEnergy;
303 double rec_NeutralEnergy = rec_NeutralHadEnergy + rec_NeutralEmEnergy;
306 std::vector <unsigned int> chMult(9, static_cast<unsigned int>(0));
307 for (
unsigned ic = 0; ic < constituents.size (); ++ic) {
308 if ( constituents[ic]->particleId() > 3 )
continue;
310 if ( trackRef.
isNull() ) {
316 unsigned int iter = 0;
317 switch (trackRef->algo()) {
319 case TrackBase::iter0:
322 case TrackBase::iter1:
325 case TrackBase::iter2:
328 case TrackBase::iter3:
331 case TrackBase::iter4:
334 case TrackBase::iter5:
337 case TrackBase::iter6:
340 case TrackBase::iter8:
347 std::cout <<
"Warning in entry " <<
entry_ <<
" : iter = " << trackRef->algo() << std::endl;
348 std::cout << ic <<
" " << *(constituents[ic]) << std::endl;
361 double cut1 = 0.0001;
362 double cut2 = 0.0001;
363 double cut3 = 0.0001;
364 double cut4 = 0.0001;
365 double cut5 = 0.0001;
366 double cut6 = 0.0001;
367 double cut7 = 0.0001;
369 double resChargedHadEnergy= 0.;
370 double resNeutralHadEnergy= 0.;
371 double resNeutralEmEnergy= 0.;
372 double resNeutralEnergy= 0.;
374 double resHCALEnergy = 0.;
375 double resNEUTEnergy = 0.;
376 if ( rec_NeutralHadEnergy > cut6 && rec_ChargedHadEnergy < cut1 ) {
377 double true_NEUTEnergy = true_NeutralHadEnergy + true_NeutralEmEnergy;
378 double true_HCALEnergy = true_NEUTEnergy - rec_NeutralEmEnergy;
379 double rec_NEUTEnergy = rec_NeutralHadEnergy+rec_NeutralEmEnergy;
380 double rec_HCALEnergy = rec_NeutralHadEnergy;
381 resHCALEnergy = (rec_HCALEnergy-true_HCALEnergy)/rec_HCALEnergy;
382 resNEUTEnergy = (rec_NEUTEnergy-true_NEUTEnergy)/rec_NEUTEnergy;
383 if ( rec_NeutralEmEnergy > cut7 ) {
391 if (true_pt > 0.0001){
392 resPt = (rec_pt -true_pt)/true_pt ;
394 if (true_ChargedHadEnergy > cut1){
395 resChargedHadEnergy = (rec_ChargedHadEnergy- true_ChargedHadEnergy)/true_ChargedHadEnergy;
397 if (true_NeutralHadEnergy > cut2){
398 resNeutralHadEnergy = (rec_NeutralHadEnergy- true_NeutralHadEnergy)/true_NeutralHadEnergy;
401 if (rec_NeutralHadEnergy > cut3){
402 resNeutralHadEnergy = (rec_NeutralHadEnergy- true_NeutralHadEnergy)/rec_NeutralHadEnergy;
404 if (true_NeutralEmEnergy > cut4){
405 resNeutralEmEnergy = (rec_NeutralEmEnergy- true_NeutralEmEnergy)/true_NeutralEmEnergy;
407 if (true_NeutralEnergy > cut5){
408 resNeutralEnergy = (rec_NeutralEnergy- true_NeutralEnergy)/true_NeutralEnergy;
415 if ( ( resPt > 0.2 && true_pt > 100. ) ||
416 ( resPt < -0.5 && true_pt > 100. ) ) {
421 <<
" resPt = " << resPt
422 <<
" resCharged " << resChargedHadEnergy
423 <<
" resNeutralHad " << resNeutralHadEnergy
424 <<
" resNeutralEm " << resNeutralEmEnergy
425 <<
" pT (T/R) " << true_pt <<
"/" << rec_pt
426 <<
" Eta (T/R) " << truth->
eta() <<
"/" << rec_eta
427 <<
" Phi (T/R) " << truth->
phi() <<
"/" << rec_phi
433 for(
unsigned j=0;
j<
pfJets.size();
j++) {
445 for(
unsigned j=0;
j<genJets.size();
j++) {
447 if ( gjo != truth &&
algo_->
deltaR(truth,gjo) < dRgo && gjo->
pt() > 0.25*truth->
pt() ) {
454 std::cout <<
"Excess probably due to overlapping jets (DR = " <<
algo_->
deltaR(genoj,pfoj) <<
"),"
455 <<
" at DeltaR(T/R) = " << dRgo <<
"/" << dRo
456 <<
" with pT(T/R) " << genoj->
pt() <<
"/" << pfoj->
pt()
457 <<
" and Eta (T/R) " << genoj->
eta() <<
"/" << pfoj->
eta()
458 <<
" and Phi (T/R) " << genoj->
phi() <<
"/" << pfoj->
phi()
467 cout <<
i <<
" =========PFJet Pt "<< rec_pt
468 <<
" eta " << rec_eta
469 <<
" phi " << rec_phi
470 <<
" Charged Had Energy " << rec_ChargedHadEnergy
471 <<
" Neutral Had Energy " << rec_NeutralHadEnergy
472 <<
" Neutral elm Energy " << rec_NeutralEmEnergy << endl;
473 cout <<
" matching Gen Jet Pt " << true_pt
474 <<
" eta " << truth->
eta()
475 <<
" phi " << truth->
phi()
476 <<
" Charged Had Energy " << true_ChargedHadEnergy
477 <<
" Neutral Had Energy " << true_NeutralHadEnergy
478 <<
" Neutral elm Energy " << true_NeutralEmEnergy << endl;
484 cout <<
"==============deltaR " << deltaR <<
" resPt " << resPt
485 <<
" resChargedHadEnergy " << resChargedHadEnergy
486 <<
" resNeutralHadEnergy " << resNeutralHadEnergy
487 <<
" resNeutralEmEnergy " << resNeutralEmEnergy
500 hNCHvsEta->Fill(true_eta, rec_ChargedMultiplicity);
510 if(plot2)
hRCHEvsEta->Fill(true_eta, resChargedHadEnergy);
511 if(plot5)
hRNeutvsEta->Fill(true_eta, resNeutralEnergy);
526 if ( pt_denom > 20. && pt_denom < 40. )
hBRPt20_40->Fill (resPt);
527 if ( pt_denom > 40. && pt_denom < 60. )
hBRPt40_60->Fill (resPt);
528 if ( pt_denom > 60. && pt_denom < 80. )
hBRPt60_80->Fill (resPt);
529 if ( pt_denom > 80. && pt_denom < 100. )
hBRPt80_100->Fill (resPt);
530 if ( pt_denom > 100. && pt_denom < 150. )
hBRPt100_150->Fill (resPt);
531 if ( pt_denom > 150. && pt_denom < 200. )
hBRPt150_200->Fill (resPt);
532 if ( pt_denom > 200. && pt_denom < 250. )
hBRPt200_250->Fill (resPt);
533 if ( pt_denom > 250. && pt_denom < 300. )
hBRPt250_300->Fill (resPt);
534 if ( pt_denom > 300. && pt_denom < 400. )
hBRPt300_400->Fill (resPt);
535 if ( pt_denom > 400. && pt_denom < 500. )
hBRPt400_500->Fill (resPt);
536 if ( pt_denom > 500. && pt_denom < 750. )
hBRPt500_750->Fill (resPt);
537 if ( pt_denom > 750. && pt_denom < 1250. )
hBRPt750_1250->Fill (resPt);
538 if ( pt_denom > 1250. && pt_denom < 2000. )
hBRPt1250_2000->Fill (resPt);
539 if ( pt_denom > 2000. && pt_denom < 5000. )
hBRPt2000_5000->Fill (resPt);
540 hBNCH->Fill(rec_ChargedMultiplicity);
549 hBNCHvsPt->Fill(pt_denom,rec_ChargedMultiplicity);
556 if(plot2)
hBRCHE->Fill(resChargedHadEnergy);
557 if(plot3)
hBRNHE->Fill(resNeutralHadEnergy);
558 if(plot4)
hBRNEE->Fill(resNeutralEmEnergy);
559 if(plot5)
hBRneut->Fill(resNeutralEnergy);
560 if(plot1)
hBRPtvsPt->Fill(pt_denom, resPt);
561 if(plot2)
hBRCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
562 if(plot3)
hBRNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
563 if(plot4)
hBRNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
564 if(plot5)
hBRneutvsPt->Fill(pt_denom, resNeutralEnergy);
579 if ( pt_denom > 20. && pt_denom < 40. )
hERPt20_40->Fill (resPt);
580 if ( pt_denom > 40. && pt_denom < 60. )
hERPt40_60->Fill (resPt);
581 if ( pt_denom > 60. && pt_denom < 80. )
hERPt60_80->Fill (resPt);
582 if ( pt_denom > 80. && pt_denom < 100. )
hERPt80_100->Fill (resPt);
583 if ( pt_denom > 100. && pt_denom < 150. )
hERPt100_150->Fill (resPt);
584 if ( pt_denom > 150. && pt_denom < 200. )
hERPt150_200->Fill (resPt);
585 if ( pt_denom > 200. && pt_denom < 250. )
hERPt200_250->Fill (resPt);
586 if ( pt_denom > 250. && pt_denom < 300. )
hERPt250_300->Fill (resPt);
587 if ( pt_denom > 300. && pt_denom < 400. )
hERPt300_400->Fill (resPt);
588 if ( pt_denom > 400. && pt_denom < 500. )
hERPt400_500->Fill (resPt);
589 if ( pt_denom > 500. && pt_denom < 750. )
hERPt500_750->Fill (resPt);
590 if ( pt_denom > 750. && pt_denom < 1250. )
hERPt750_1250->Fill (resPt);
591 if ( pt_denom > 1250. && pt_denom < 2000. )
hERPt1250_2000->Fill (resPt);
592 if ( pt_denom > 2000. && pt_denom < 5000. )
hERPt2000_5000->Fill (resPt);
593 hENCH->Fill(rec_ChargedMultiplicity);
594 hENCHvsPt->Fill(pt_denom,rec_ChargedMultiplicity);
609 if(plot2)
hERCHE->Fill(resChargedHadEnergy);
610 if(plot3)
hERNHE->Fill(resNeutralHadEnergy);
611 if(plot4)
hERNEE->Fill(resNeutralEmEnergy);
612 if(plot5)
hERneut->Fill(resNeutralEnergy);
613 if(plot1)
hERPtvsPt->Fill(pt_denom, resPt);
614 if(plot2)
hERCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
615 if(plot3)
hERNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
616 if(plot4)
hERNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
617 if(plot5)
hERneutvsPt->Fill(pt_denom, resNeutralEnergy);
631 if ( pt_denom > 20. && pt_denom < 40. )
hFRPt20_40->Fill (resPt);
632 if ( pt_denom > 40. && pt_denom < 60. )
hFRPt40_60->Fill (resPt);
633 if ( pt_denom > 60. && pt_denom < 80. )
hFRPt60_80->Fill (resPt);
634 if ( pt_denom > 80. && pt_denom < 100. )
hFRPt80_100->Fill (resPt);
635 if ( pt_denom > 100. && pt_denom < 150. )
hFRPt100_150->Fill (resPt);
636 if ( pt_denom > 150. && pt_denom < 200. )
hFRPt150_200->Fill (resPt);
637 if ( pt_denom > 200. && pt_denom < 250. )
hFRPt200_250->Fill (resPt);
638 if ( pt_denom > 250. && pt_denom < 300. )
hFRPt250_300->Fill (resPt);
639 if ( pt_denom > 300. && pt_denom < 400. )
hFRPt300_400->Fill (resPt);
640 if ( pt_denom > 400. && pt_denom < 500. )
hFRPt400_500->Fill (resPt);
641 if ( pt_denom > 500. && pt_denom < 750. )
hFRPt500_750->Fill (resPt);
642 if ( pt_denom > 750. && pt_denom < 1250. )
hFRPt750_1250->Fill (resPt);
643 if ( pt_denom > 1250. && pt_denom < 2000. )
hFRPt1250_2000->Fill (resPt);
644 if ( pt_denom > 2000. && pt_denom < 5000. )
hFRPt2000_5000->Fill (resPt);
651 if(plot2)
hFRCHE->Fill(resChargedHadEnergy);
652 if(plot3)
hFRNHE->Fill(resNeutralHadEnergy);
653 if(plot4)
hFRNEE->Fill(resNeutralEmEnergy);
654 if(plot5)
hFRneut->Fill(resNeutralEnergy);
655 if(plot1)
hFRPtvsPt->Fill(pt_denom, resPt);
656 if(plot2)
hFRCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
657 if(plot3)
hFRNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
658 if(plot4)
hFRNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
659 if(plot5)
hFRneutvsPt->Fill(pt_denom, resNeutralEnergy);
virtual double p() const
magnitude of momentum vector
virtual float pt() const
transverse momentum
virtual float phi() const
momentum azimuthal angle
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.
bool isNull() const
Checks for null.
virtual float eta() const
momentum pseudorapidity
static double deltaR(const T *, const U *)
Abs< T >::type abs(const T &t)
Jets made from MC generator particles.
double deltaR(double eta1, double eta2, double phi1, double phi2)
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_
void printGenJet(const reco::GenJet *)
float chargedHadronEnergy() const
chargedHadronEnergy