236 for(
unsigned i=0;
i<
pfJets.size();
i++) {
239 unsigned highJets = 0;
240 for(
unsigned j=0;
j<
pfJets.size();
j++) {
247 double rec_pt = pfj.
pt();
248 double rec_eta = pfj.
eta();
249 double rec_phi = pfj.
phi();
266 bool Forward =
false;
267 if (
std::abs(rec_eta) < 1.4 ) Barrel =
true;
285 double true_E = truth->
p();
286 double true_pt = truth->
pt();
287 double true_eta = truth->
eta();
288 double true_phi = truth->
phi();
291 else {pt_denom = true_pt;}
293 double true_ChargedHadEnergy;
294 double true_NeutralHadEnergy;
295 double true_NeutralEmEnergy;
296 gettrue (truth, true_ChargedHadEnergy, true_NeutralHadEnergy, true_NeutralEmEnergy);
297 double true_NeutralEnergy = true_NeutralHadEnergy + true_NeutralEmEnergy;
301 double rec_NeutralEnergy = rec_NeutralHadEnergy + rec_NeutralEmEnergy;
304 std::vector <unsigned int> chMult(7, static_cast<unsigned int>(0));
305 for (
unsigned ic = 0; ic < constituents.size (); ++ic) {
306 if ( constituents[ic]->particleId() > 3 )
continue;
308 if ( trackRef.
isNull() ) {
314 unsigned int iter = 0;
315 switch (trackRef->algo()) {
317 case TrackBase::iter0:
320 case TrackBase::iter1:
323 case TrackBase::iter2:
326 case TrackBase::iter3:
329 case TrackBase::iter4:
332 case TrackBase::iter5:
337 std::cout <<
"Warning in entry " <<
entry_ <<
" : iter = 6... " << std::endl;
350 double cut1 = 0.0001;
351 double cut2 = 0.0001;
352 double cut3 = 0.0001;
353 double cut4 = 0.0001;
354 double cut5 = 0.0001;
355 double cut6 = 0.0001;
356 double cut7 = 0.0001;
358 double resChargedHadEnergy= 0.;
359 double resNeutralHadEnergy= 0.;
360 double resNeutralEmEnergy= 0.;
361 double resNeutralEnergy= 0.;
363 double resHCALEnergy = 0.;
364 double resNEUTEnergy = 0.;
365 if ( rec_NeutralHadEnergy > cut6 && rec_ChargedHadEnergy < cut1 ) {
366 double true_NEUTEnergy = true_NeutralHadEnergy + true_NeutralEmEnergy;
367 double true_HCALEnergy = true_NEUTEnergy - rec_NeutralEmEnergy;
368 double rec_NEUTEnergy = rec_NeutralHadEnergy+rec_NeutralEmEnergy;
369 double rec_HCALEnergy = rec_NeutralHadEnergy;
370 resHCALEnergy = (rec_HCALEnergy-true_HCALEnergy)/rec_HCALEnergy;
371 resNEUTEnergy = (rec_NEUTEnergy-true_NEUTEnergy)/rec_NEUTEnergy;
372 if ( rec_NeutralEmEnergy > cut7 ) {
380 if (true_pt > 0.0001){
381 resPt = (rec_pt -true_pt)/true_pt ;
383 if (true_ChargedHadEnergy > cut1){
384 resChargedHadEnergy = (rec_ChargedHadEnergy- true_ChargedHadEnergy)/true_ChargedHadEnergy;
386 if (true_NeutralHadEnergy > cut2){
387 resNeutralHadEnergy = (rec_NeutralHadEnergy- true_NeutralHadEnergy)/true_NeutralHadEnergy;
390 if (rec_NeutralHadEnergy > cut3){
391 resNeutralHadEnergy = (rec_NeutralHadEnergy- true_NeutralHadEnergy)/rec_NeutralHadEnergy;
393 if (true_NeutralEmEnergy > cut4){
394 resNeutralEmEnergy = (rec_NeutralEmEnergy- true_NeutralEmEnergy)/true_NeutralEmEnergy;
396 if (true_NeutralEnergy > cut5){
397 resNeutralEnergy = (rec_NeutralEnergy- true_NeutralEnergy)/true_NeutralEnergy;
404 if ( ( resPt > 0.2 && true_pt > 100. ) ||
405 ( resPt < -0.5 && true_pt > 100. ) ) {
410 <<
" resPt = " << resPt
411 <<
" resCharged " << resChargedHadEnergy
412 <<
" resNeutralHad " << resNeutralHadEnergy
413 <<
" resNeutralEm " << resNeutralEmEnergy
414 <<
" pT (T/R) " << true_pt <<
"/" << rec_pt
415 <<
" Eta (T/R) " << truth->
eta() <<
"/" << rec_eta
416 <<
" Phi (T/R) " << truth->
phi() <<
"/" << rec_phi
422 for(
unsigned j=0;
j<
pfJets.size();
j++) {
434 for(
unsigned j=0;
j<genJets.size();
j++) {
436 if ( gjo != truth &&
algo_->
deltaR(truth,gjo) < dRgo && gjo->
pt() > 0.25*truth->
pt() ) {
443 std::cout <<
"Excess probably due to overlapping jets (DR = " <<
algo_->
deltaR(genoj,pfoj) <<
"),"
444 <<
" at DeltaR(T/R) = " << dRgo <<
"/" << dRo
445 <<
" with pT(T/R) " << genoj->
pt() <<
"/" << pfoj->
pt()
446 <<
" and Eta (T/R) " << genoj->
eta() <<
"/" << pfoj->
eta()
447 <<
" and Phi (T/R) " << genoj->
phi() <<
"/" << pfoj->
phi()
456 cout <<
i <<
" =========PFJet Pt "<< rec_pt
457 <<
" eta " << rec_eta
458 <<
" phi " << rec_phi
459 <<
" Charged Had Energy " << rec_ChargedHadEnergy
460 <<
" Neutral Had Energy " << rec_NeutralHadEnergy
461 <<
" Neutral elm Energy " << rec_NeutralEmEnergy << endl;
462 cout <<
" matching Gen Jet Pt " << true_pt
463 <<
" eta " << truth->
eta()
464 <<
" phi " << truth->
phi()
465 <<
" Charged Had Energy " << true_ChargedHadEnergy
466 <<
" Neutral Had Energy " << true_NeutralHadEnergy
467 <<
" Neutral elm Energy " << true_NeutralEmEnergy << endl;
473 cout <<
"==============deltaR " << deltaR <<
" resPt " << resPt
474 <<
" resChargedHadEnergy " << resChargedHadEnergy
475 <<
" resNeutralHadEnergy " << resNeutralHadEnergy
476 <<
" resNeutralEmEnergy " << resNeutralEmEnergy
489 hNCHvsEta->Fill(true_eta, rec_ChargedMultiplicity);
498 if(plot2)
hRCHEvsEta->Fill(true_eta, resChargedHadEnergy);
499 if(plot5)
hRNeutvsEta->Fill(true_eta, resNeutralEnergy);
514 if ( pt_denom > 20. && pt_denom < 40. )
hBRPt20_40->Fill (resPt);
515 if ( pt_denom > 40. && pt_denom < 60. )
hBRPt40_60->Fill (resPt);
516 if ( pt_denom > 60. && pt_denom < 80. )
hBRPt60_80->Fill (resPt);
517 if ( pt_denom > 80. && pt_denom < 100. )
hBRPt80_100->Fill (resPt);
518 if ( pt_denom > 100. && pt_denom < 150. )
hBRPt100_150->Fill (resPt);
519 if ( pt_denom > 150. && pt_denom < 200. )
hBRPt150_200->Fill (resPt);
520 if ( pt_denom > 200. && pt_denom < 250. )
hBRPt200_250->Fill (resPt);
521 if ( pt_denom > 250. && pt_denom < 300. )
hBRPt250_300->Fill (resPt);
522 if ( pt_denom > 300. && pt_denom < 400. )
hBRPt300_400->Fill (resPt);
523 if ( pt_denom > 400. && pt_denom < 500. )
hBRPt400_500->Fill (resPt);
524 if ( pt_denom > 500. && pt_denom < 750. )
hBRPt500_750->Fill (resPt);
525 if ( pt_denom > 750. && pt_denom < 1250. )
hBRPt750_1250->Fill (resPt);
526 if ( pt_denom > 1250. && pt_denom < 2000. )
hBRPt1250_2000->Fill (resPt);
527 if ( pt_denom > 2000. && pt_denom < 5000. )
hBRPt2000_5000->Fill (resPt);
528 hBNCH->Fill(rec_ChargedMultiplicity);
536 hBNCHvsPt->Fill(pt_denom,rec_ChargedMultiplicity);
543 if(plot2)
hBRCHE->Fill(resChargedHadEnergy);
544 if(plot3)
hBRNHE->Fill(resNeutralHadEnergy);
545 if(plot4)
hBRNEE->Fill(resNeutralEmEnergy);
546 if(plot5)
hBRneut->Fill(resNeutralEnergy);
547 if(plot1)
hBRPtvsPt->Fill(pt_denom, resPt);
548 if(plot2)
hBRCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
549 if(plot3)
hBRNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
550 if(plot4)
hBRNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
551 if(plot5)
hBRneutvsPt->Fill(pt_denom, resNeutralEnergy);
566 if ( pt_denom > 20. && pt_denom < 40. )
hERPt20_40->Fill (resPt);
567 if ( pt_denom > 40. && pt_denom < 60. )
hERPt40_60->Fill (resPt);
568 if ( pt_denom > 60. && pt_denom < 80. )
hERPt60_80->Fill (resPt);
569 if ( pt_denom > 80. && pt_denom < 100. )
hERPt80_100->Fill (resPt);
570 if ( pt_denom > 100. && pt_denom < 150. )
hERPt100_150->Fill (resPt);
571 if ( pt_denom > 150. && pt_denom < 200. )
hERPt150_200->Fill (resPt);
572 if ( pt_denom > 200. && pt_denom < 250. )
hERPt200_250->Fill (resPt);
573 if ( pt_denom > 250. && pt_denom < 300. )
hERPt250_300->Fill (resPt);
574 if ( pt_denom > 300. && pt_denom < 400. )
hERPt300_400->Fill (resPt);
575 if ( pt_denom > 400. && pt_denom < 500. )
hERPt400_500->Fill (resPt);
576 if ( pt_denom > 500. && pt_denom < 750. )
hERPt500_750->Fill (resPt);
577 if ( pt_denom > 750. && pt_denom < 1250. )
hERPt750_1250->Fill (resPt);
578 if ( pt_denom > 1250. && pt_denom < 2000. )
hERPt1250_2000->Fill (resPt);
579 if ( pt_denom > 2000. && pt_denom < 5000. )
hERPt2000_5000->Fill (resPt);
580 hENCH->Fill(rec_ChargedMultiplicity);
581 hENCHvsPt->Fill(pt_denom,rec_ChargedMultiplicity);
595 if(plot2)
hERCHE->Fill(resChargedHadEnergy);
596 if(plot3)
hERNHE->Fill(resNeutralHadEnergy);
597 if(plot4)
hERNEE->Fill(resNeutralEmEnergy);
598 if(plot5)
hERneut->Fill(resNeutralEnergy);
599 if(plot1)
hERPtvsPt->Fill(pt_denom, resPt);
600 if(plot2)
hERCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
601 if(plot3)
hERNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
602 if(plot4)
hERNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
603 if(plot5)
hERneutvsPt->Fill(pt_denom, resNeutralEnergy);
617 if ( pt_denom > 20. && pt_denom < 40. )
hFRPt20_40->Fill (resPt);
618 if ( pt_denom > 40. && pt_denom < 60. )
hFRPt40_60->Fill (resPt);
619 if ( pt_denom > 60. && pt_denom < 80. )
hFRPt60_80->Fill (resPt);
620 if ( pt_denom > 80. && pt_denom < 100. )
hFRPt80_100->Fill (resPt);
621 if ( pt_denom > 100. && pt_denom < 150. )
hFRPt100_150->Fill (resPt);
622 if ( pt_denom > 150. && pt_denom < 200. )
hFRPt150_200->Fill (resPt);
623 if ( pt_denom > 200. && pt_denom < 250. )
hFRPt200_250->Fill (resPt);
624 if ( pt_denom > 250. && pt_denom < 300. )
hFRPt250_300->Fill (resPt);
625 if ( pt_denom > 300. && pt_denom < 400. )
hFRPt300_400->Fill (resPt);
626 if ( pt_denom > 400. && pt_denom < 500. )
hFRPt400_500->Fill (resPt);
627 if ( pt_denom > 500. && pt_denom < 750. )
hFRPt500_750->Fill (resPt);
628 if ( pt_denom > 750. && pt_denom < 1250. )
hFRPt750_1250->Fill (resPt);
629 if ( pt_denom > 1250. && pt_denom < 2000. )
hFRPt1250_2000->Fill (resPt);
630 if ( pt_denom > 2000. && pt_denom < 5000. )
hFRPt2000_5000->Fill (resPt);
637 if(plot2)
hFRCHE->Fill(resChargedHadEnergy);
638 if(plot3)
hFRNHE->Fill(resNeutralHadEnergy);
639 if(plot4)
hFRNEE->Fill(resNeutralEmEnergy);
640 if(plot5)
hFRneut->Fill(resNeutralEnergy);
641 if(plot1)
hFRPtvsPt->Fill(pt_denom, resPt);
642 if(plot2)
hFRCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
643 if(plot3)
hFRNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
644 if(plot4)
hFRNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
645 if(plot5)
hFRneutvsPt->Fill(pt_denom, resNeutralEnergy);
virtual double p() const
magnitude of momentum vector
int chargedMultiplicity() const
chargedMultiplicity
Jets made from PFObjects.
virtual double eta() const
momentum pseudorapidity
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.
static double deltaR(const T *, const U *)
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 *)
virtual double pt() const
transverse momentum
double resNeutralHadEnergyMax_
void printPFJet(const reco::PFJet *)
virtual std::vector< reco::PFCandidatePtr > getPFConstituents() const
get all constituents
float neutralHadronEnergy() const
neutralHadronEnergy
double resNeutralEmEnergyMax_
virtual double phi() const
momentum azimuthal angle
void printGenJet(const reco::GenJet *)
float chargedHadronEnergy() const
chargedHadronEnergy