3175 using namespace edm;
3179 const float BARL = 1.4442;
3181 const float END_HI = 2.5;
3186 e.getByLabel(
"trackAssociatorByHitsForPhotonValidation", theHitsAssociator);
3190 LogInfo(
"PhotonValidator") <<
"PhotonValidator Analyzing event number: " <<
e.id() <<
" Global Counter " <<
nEvt_
3208 if (!photonHandle.
isValid()) {
3209 edm::LogError(
"PhotonProducer") <<
"Error! Can't get the Photon collection " << std::endl;
3216 if (!pfCandidateHandle.
isValid()) {
3217 edm::LogError(
"PhotonValidator") <<
"Error! Can't get the product pfCandidates " << std::endl;
3222 if (
fName_ ==
"pfPhotonValidator") {
3224 if (!phoToParticleBasedIsoMapHandle.
isValid()) {
3225 edm::LogInfo(
"PhotonValidator") <<
"Error! Can't get the product: valueMap photons to particle based iso "
3228 phoToParticleBasedIsoMap = *(phoToParticleBasedIsoMapHandle.
product());
3246 if ((*itHits)->isValid()) {
3261 if ((*itHits)->isValid()) {
3272 std::vector<SimTrack> theSimTracks;
3273 std::vector<SimVertex> theSimVertices;
3285 theSimTracks.insert(theSimTracks.end(), SimTk->begin(), SimTk->end());
3286 theSimVertices.insert(theSimVertices.end(), SimVtx->begin(), SimVtx->end());
3313 std::vector<reco::PhotonCollection::const_iterator> StoRMatchedConvertedPhotons;
3323 OISimToReco =
trackAssociator->associateSimToReco(outInTrkHandle, ElectronTPHandle);
3324 IOSimToReco =
trackAssociator->associateSimToReco(inOutTrkHandle, ElectronTPHandle);
3326 OIRecoToSim =
trackAssociator->associateRecoToSim(outInTrkHandle, ElectronTPHandle);
3327 IORecoToSim =
trackAssociator->associateRecoToSim(inOutTrkHandle, ElectronTPHandle);
3331 vector<reco::SimToRecoCollection*> StoRCollPtrs;
3332 StoRCollPtrs.push_back(&OISimToReco);
3333 StoRCollPtrs.push_back(&IOSimToReco);
3334 vector<reco::RecoToSimCollection*> RtoSCollPtrs;
3335 RtoSCollPtrs.push_back(&OIRecoToSim);
3336 RtoSCollPtrs.push_back(&IORecoToSim);
3338 for (
int i = 0;
i < 2;
i++)
3340 for (
int i = 0;
i < 2;
i++)
3343 std::vector<reco::PhotonRef> myPhotons;
3345 for (
unsigned int iPho = 0; iPho < photonHandle->size(); iPho++) {
3348 if (fabs(phoRef->eta()) > 2.5)
3350 myPhotons.push_back(phoRef);
3353 std::sort(myPhotons.begin(), myPhotons.end(), sortPhotons());
3355 if (myPhotons.size() >= 2) {
3356 if (myPhotons[0]->
et() > 40 && myPhotons[1]->et() > 25) {
3362 float gamgamMass2 = p12.Dot(p12);
3363 float gamgamMass2_regr1 = p12_regr1.Dot(p12_regr1);
3364 float gamgamMass2_regr2 = p12_regr2.Dot(p12_regr2);
3367 if (gamgamMass2 > 0) {
3370 if (myPhotons[0]->isEB() && myPhotons[1]->isEB())
3372 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3373 (myPhotons[0]->isEB() && myPhotons[1]->isEE()))
3376 if (myPhotons[0]->
r9() > 0.94 && myPhotons[1]->
r9() > 0.94) {
3378 if (myPhotons[0]->isEB() && myPhotons[1]->isEB())
3380 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3381 (myPhotons[0]->isEB() && myPhotons[1]->isEE()))
3391 if (chi2Prob1 > 0.0005 && chi2Prob2 > 0.0005) {
3393 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3396 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3397 (myPhotons[0]->isEB() && myPhotons[1]->isEE())) {
3403 myPhotons[1]->
r9() > 0.93) {
3407 if (chi2Prob1 > 0.0005) {
3409 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3412 if (myPhotons[0]->isEE() || myPhotons[1]->isEE()) {
3418 myPhotons[0]->
r9() > 0.93) {
3422 if (chi2Prob1 > 0.0005) {
3424 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3427 if (myPhotons[0]->isEE() || myPhotons[1]->isEE()) {
3436 if (gamgamMass2_regr1 > 0) {
3439 if (myPhotons[0]->isEB() && myPhotons[1]->isEB())
3441 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3442 (myPhotons[0]->isEB() && myPhotons[1]->isEE()))
3445 if (myPhotons[0]->
r9() > 0.94 && myPhotons[1]->
r9() > 0.94) {
3447 if (myPhotons[0]->isEB() && myPhotons[1]->isEB())
3449 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3450 (myPhotons[0]->isEB() && myPhotons[1]->isEE()))
3461 if (chi2Prob1 > 0.0005 && chi2Prob2 > 0.0005) {
3463 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3466 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3467 (myPhotons[0]->isEB() && myPhotons[1]->isEE())) {
3473 myPhotons[1]->
r9() > 0.93) {
3477 if (chi2Prob1 > 0.0005) {
3479 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3482 if (myPhotons[0]->isEE() || myPhotons[1]->isEE()) {
3488 myPhotons[0]->
r9() > 0.93) {
3492 if (chi2Prob1 > 0.0005) {
3494 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3497 if (myPhotons[0]->isEE() || myPhotons[1]->isEE()) {
3506 if (gamgamMass2_regr2 > 0) {
3509 if (myPhotons[0]->isEB() && myPhotons[1]->isEB())
3511 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3512 (myPhotons[0]->isEB() && myPhotons[1]->isEE()))
3515 if (myPhotons[0]->
r9() > 0.94 && myPhotons[1]->
r9() > 0.94) {
3517 if (myPhotons[0]->isEB() && myPhotons[1]->isEB())
3519 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3520 (myPhotons[0]->isEB() && myPhotons[1]->isEE()))
3531 if (chi2Prob1 > 0.0005 && chi2Prob2 > 0.0005) {
3533 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3536 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3537 (myPhotons[0]->isEB() && myPhotons[1]->isEE())) {
3543 myPhotons[1]->
r9() > 0.93) {
3547 if (chi2Prob1 > 0.0005) {
3549 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3552 if (myPhotons[0]->isEE() || myPhotons[1]->isEE()) {
3558 myPhotons[0]->
r9() > 0.93) {
3562 if (chi2Prob1 > 0.0005) {
3564 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3567 if (myPhotons[0]->isEE() || myPhotons[1]->isEE()) {
3578 for (std::vector<PhotonMCTruth>::const_iterator mcPho = mcPhotons.begin(); mcPho != mcPhotons.end(); mcPho++) {
3582 for (HepMC::GenEvent::particle_const_iterator mcIter = myGenEvent->particles_begin();
3583 mcIter != myGenEvent->particles_end();
3585 if ((*mcIter)->pdg_id() != 22)
3587 bool isTheSame =
false;
3589 if ((*mcIter)->production_vertex()) {
3590 if ((*mcIter)->production_vertex()->particles_begin(
HepMC::parents) !=
3592 mother = *((*mcIter)->production_vertex()->particles_begin(
HepMC::parents));
3595 float mcPhi = (*mcPho).fourMomentum().phi();
3597 mcEta_ = (*mcPho).fourMomentum().pseudoRapidity();
3599 mcConvR_ = (*mcPho).vertex().perp();
3606 if (fabs(
mcEta_) > END_HI)
3609 if (mother ==
nullptr || (mother !=
nullptr && mother->pdg_id() == 22) ||
3610 (mother !=
nullptr && mother->pdg_id() == 25) || (mother !=
nullptr && mother->pdg_id() == 35)) {
3611 double dPt = fabs((*mcIter)->momentum().perp() - (*mcPho).fourMomentum().et());
3612 float phiMother = (*mcIter)->momentum().phi();
3614 double dEta = fabs((*mcIter)->momentum().eta() - (*mcPho).fourMomentum().pseudoRapidity());
3616 if (
dEta <= 0.0001 &&
dPhi <= 0.0001 && dPt <= 0.0001)
3634 bool goodSimConversion =
false;
3635 bool visibleConversion =
false;
3636 bool visibleConversionsWithTwoSimTracks =
false;
3637 if ((*mcPho).isAConversion() == 1) {
3651 (fabs(
mcEta_) > BARL && fabs(
mcEta_) <= END_HI && fabs((*mcPho).vertex().z()) < 210))
3652 visibleConversion =
true;
3657 if (fabs(
tp->vx() - (*mcPho).vertex().x()) < 0.001 && fabs(
tp->vy() - (*mcPho).vertex().y()) < 0.001 &&
3658 fabs(
tp->vz() - (*mcPho).vertex().z()) < 0.001) {
3663 visibleConversionsWithTwoSimTracks =
true;
3664 goodSimConversion =
false;
3666 if (visibleConversion && visibleConversionsWithTwoSimTracks)
3667 goodSimConversion =
true;
3668 if (goodSimConversion) {
3689 float minDelta = 10000.;
3690 std::vector<reco::PhotonRef> thePhotons;
3695 for (
unsigned int iPho = 0; iPho < photonHandle->size(); iPho++) {
3697 thePhotons.push_back(aPho);
3698 float phiPho = aPho->phi();
3699 float etaPho = aPho->eta();
3732 bool phoIsInBarrel =
false;
3733 bool phoIsInEndcap =
false;
3734 bool phoIsInEndcapP =
false;
3735 bool phoIsInEndcapM =
false;
3739 if (fabs(matchingPho->superCluster()->position().eta()) < 1.479) {
3740 phoIsInBarrel =
true;
3742 phoIsInEndcap =
true;
3743 if (matchingPho->superCluster()->position().eta() > 0)
3744 phoIsInEndcapP =
true;
3745 if (matchingPho->superCluster()->position().eta() < 0)
3746 phoIsInEndcapM =
true;
3750 if (phoIsInBarrel) {
3753 if (!ecalRecHitHandle.
isValid()) {
3756 edm::LogError(
"PhotonProducer") <<
"Error! Can't get the product " <<
l.module;
3760 }
else if (phoIsInEndcap) {
3763 if (!ecalRecHitHandle.
isValid()) {
3766 edm::LogError(
"PhotonProducer") <<
"Error! Can't get the product " <<
l.module;
3773 float photonE = matchingPho->energy();
3774 float sigmaEoE = matchingPho->getCorrectedEnergyError(matchingPho->getCandidateP4type()) / matchingPho->energy();
3776 float photonEt = matchingPho->pt();
3779 float r9 = matchingPho->r9();
3781 float r1 = matchingPho->r1x5();
3782 float r2 = matchingPho->r2x5();
3783 float sigmaIetaIeta = matchingPho->sigmaIetaIeta();
3785 float hOverE = matchingPho->hadronicOverEm();
3786 float newhOverE = matchingPho->hadTowOverEm();
3787 float ecalIso = matchingPho->ecalRecHitSumEtConeDR04();
3788 float hcalIso = matchingPho->hcalTowerSumEtConeDR04();
3789 float newhcalIso = matchingPho->hcalTowerSumEtBcConeDR04();
3790 float trkIso = matchingPho->trkSumPtSolidConeDR04();
3791 float nIsoTrk = matchingPho->nTrkSolidConeDR04();
3793 float chargedHadIso = matchingPho->chargedHadronIso();
3794 float neutralHadIso = matchingPho->neutralHadronIso();
3795 float photonIso = matchingPho->photonIso();
3796 float etOutsideMustache = matchingPho->etOutsideMustache();
3797 int nClusterOutsideMustache = matchingPho->nClusterOutsideMustache();
3798 float pfMVA = matchingPho->pfMVA();
3800 std::vector<std::pair<DetId, float> >::const_iterator rhIt;
3801 bool atLeastOneDeadChannel =
false;
3803 bcIt != matchingPho->superCluster()->clustersEnd();
3805 for (rhIt = (*bcIt)->hitsAndFractions().begin(); rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) {
3808 if (rhIt->first == (*it).id()) {
3809 if ((*it).recoFlag() == 9) {
3810 atLeastOneDeadChannel =
true;
3818 if (atLeastOneDeadChannel) {
3836 h_scEt_[
type][0]->
Fill(matchingPho->superCluster()->energy() / cosh(matchingPho->superCluster()->eta()));
3838 h_psE_->
Fill(matchingPho->superCluster()->preshowerEnergy());
3911 h_phoDEta_[0]->
Fill(matchingPho->eta() - (*mcPho).fourMomentum().eta());
3915 h_nConv_[0][0]->
Fill(
float(matchingPho->conversions().size()));
3916 h_nConv_[1][0]->
Fill(
float(matchingPho->conversionsOneLeg().size()));
3926 p_eResVsEt_[0][0]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
3931 h2_eResVsEt_[0][0]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
3935 h2_sceResVsR9_[0]->
Fill(
r9, matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
3939 p_sceResVsR9_[0]->
Fill(
r9, matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
3941 if ((*mcPho).isAConversion() == 0) {
3972 if (photonE / (*mcPho).fourMomentum().e() < 0.3 && photonE / (*mcPho).fourMomentum().e() > 0.1) {
3975 if ((
r9 > 0.94 && phoIsInBarrel) || (
r9 > 0.95 && phoIsInEndcap)) {
3981 h2_eResVsEt_[0][1]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
3982 p_eResVsEt_[0][1]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
3987 }
else if ((
r9 <= 0.94 && phoIsInBarrel) || (
r9 <= 0.95 && phoIsInEndcap)) {
3992 p_eResVsEt_[0][2]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
3998 h2_eResVsEt_[0][2]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4003 if (phoIsInBarrel) {
4005 h_scEt_[
type][1]->
Fill(matchingPho->superCluster()->energy() / cosh(matchingPho->superCluster()->eta()));
4029 h_nConv_[1][1]->
Fill(
float(matchingPho->conversionsOneLeg().size()));
4035 p_sceResVsR9_[1]->
Fill(
r9, matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4038 h2_sceResVsR9_[1]->
Fill(
r9, matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4041 h2_eResVsEt_[1][0]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4043 p_eResVsEt_[1][0]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4044 p_eResVsNVtx_[1][0]->
Fill(
float(vtxH->size()), photonE / (*mcPho).fourMomentum().e());
4054 h2_eResVsEt_[1][1]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4055 p_eResVsEt_[1][1]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4056 p_eResVsNVtx_[1][1]->
Fill(
float(vtxH->size()), photonE / (*mcPho).fourMomentum().e());
4065 p_eResVsEt_[1][2]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4066 p_eResVsNVtx_[1][2]->
Fill(
float(vtxH->size()), photonE / (*mcPho).fourMomentum().e());
4070 h2_eResVsEt_[1][2]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4075 if (phoIsInEndcap) {
4077 h_scEt_[
type][2]->
Fill(matchingPho->superCluster()->energy() / cosh(matchingPho->superCluster()->eta()));
4101 h_nConv_[1][2]->
Fill(
float(matchingPho->conversionsOneLeg().size()));
4107 p_sceResVsR9_[2]->
Fill(
r9, matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4110 h2_sceResVsR9_[2]->
Fill(
r9, matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4113 h2_eResVsEt_[2][0]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4116 p_eResVsEt_[2][0]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4117 p_eResVsNVtx_[2][0]->
Fill(
float(vtxH->size()), photonE / (*mcPho).fourMomentum().e());
4127 h2_eResVsEt_[2][1]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4128 p_eResVsEt_[2][1]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4129 p_eResVsNVtx_[2][1]->
Fill(
float(vtxH->size()), photonE / (*mcPho).fourMomentum().e());
4138 p_eResVsEt_[2][2]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4139 p_eResVsNVtx_[2][2]->
Fill(
float(vtxH->size()), photonE / (*mcPho).fourMomentum().e());
4144 h2_eResVsEt_[2][2]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4151 if (
fName_ ==
"pfPhotonValidator") {
4152 float SumPtIsoValCh = 0.;
4153 float SumPtIsoValNh = 0.;
4154 float SumPtIsoValPh = 0.;
4156 float SumPtIsoValCleanCh = 0.;
4157 float SumPtIsoValCleanNh = 0.;
4158 float SumPtIsoValCleanPh = 0.;
4160 for (
unsigned int lCand = 0; lCand < pfCandidateHandle->size(); lCand++) {
4162 float dR =
deltaR(matchingPho->eta(), matchingPho->phi(), pfCandRef->eta(), pfCandRef->phi());
4172 SumPtIsoValCh += pfCandRef->pt();
4180 SumPtIsoValNh += pfCandRef->pt();
4188 SumPtIsoValPh += pfCandRef->pt();
4197 for (std::vector<reco::PFCandidateRef>::const_iterator
i = phoToParticleBasedIsoMap[matchingPho].begin();
4198 i != phoToParticleBasedIsoMap[matchingPho].
end();
4200 if ((*
i) == pfCandRef) {
4208 SumPtIsoValCleanCh += pfCandRef->pt();
4216 SumPtIsoValCleanNh += pfCandRef->pt();
4224 SumPtIsoValCleanPh += pfCandRef->pt();
4241 if (phoIsInBarrel) {
4259 if (!(visibleConversion && visibleConversionsWithTwoSimTracks))
4270 if (fabs(
mcEta_) <= 1.) {
4280 bool atLeastOneRecoTwoTrackConversion =
false;
4281 for (
unsigned int iConv = 0; iConv <
conversions.size(); iConv++) {
4283 double like = aConv->MVAout();
4292 const std::vector<edm::RefToBase<reco::Track> >
tracks = aConv->tracks();
4295 atLeastOneRecoTwoTrackConversion =
true;
4299 if (
tracks.size() == 2) {
4314 std::map<const reco::Track*, TrackingParticleRef> myAss;
4315 std::map<const reco::Track*, TrackingParticleRef>::const_iterator itAss;
4316 std::map<reco::TrackRef, TrackingParticleRef>::const_iterator itAssMin;
4317 std::map<reco::TrackRef, TrackingParticleRef>::const_iterator itAssMax;
4321 for (
unsigned int i = 0;
i <
tracks.size();
i++) {
4338 std::vector<std::pair<RefToBase<reco::Track>,
double> > trackV;
4342 trackV = (std::vector<std::pair<RefToBase<reco::Track>,
double> >)
q[
theConvTP_[0]];
4344 trackV = (std::vector<std::pair<RefToBase<reco::Track>,
double> >)
q[
theConvTP_[1]];
4358 float refP = -99999.;
4359 float refPt = -99999.;
4360 if (aConv->conversionVertex().isValid()) {
4361 refP =
sqrt(aConv->refittedPairMomentum().Mag2());
4362 refPt =
sqrt(aConv->refittedPairMomentum().perp2());
4364 float invM = aConv->pairInvariantMass();
4373 if (
tracks.size() == 2) {
4380 if (!aConv->caloCluster().empty())
4383 float trkProvenance = 3;
4384 if (
tracks[0]->
algoName() ==
"outInEcalSeededConv" &&
tracks[1]->algoName() ==
"outInEcalSeededConv")
4386 if (
tracks[0]->
algoName() ==
"inOutEcalSeededConv" &&
tracks[1]->algoName() ==
"inOutEcalSeededConv")
4389 (
tracks[1]->algoName() ==
"outInEcalSeededConv" &&
tracks[0]->algoName() ==
"inOutEcalSeededConv"))
4391 if (trkProvenance == 3) {
4408 if (!aConv->caloCluster().empty())
4418 if (aConv->conversionVertex().isValid()) {
4428 if (chi2Prob > 0.0005) {
4436 if (chi2Prob > 0.0005) {
4437 if (!aConv->caloCluster().empty()) {
4440 h_convERes_[0][0]->
Fill(aConv->caloCluster()[0]->energy() / (*mcPho).fourMomentum().e());
4450 if (phoIsInBarrel) {
4451 if (!aConv->caloCluster().empty())
4452 h_convERes_[0][1]->
Fill(aConv->caloCluster()[0]->energy() / (*mcPho).fourMomentum().e());
4457 if (phoIsInEndcap) {
4458 if (!aConv->caloCluster().empty())
4459 h_convERes_[0][2]->
Fill(aConv->caloCluster()[0]->energy() / (*mcPho).fourMomentum().e());
4473 float eoverp = -99999.;
4475 if (aConv->conversionVertex().isValid()) {
4476 eoverp = photonE /
sqrt(aConv->refittedPairMomentum().Mag2());
4483 matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4500 matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4503 matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4511 float dPhiTracksAtVtx = aConv->dPhiTracksAtVtx();
4528 if (phoIsInBarrel) {
4530 if (aConv->conversionVertex().isValid()) {
4542 eoverp, matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4550 if (phoIsInEndcap) {
4552 if (aConv->conversionVertex().isValid()) {
4558 eoverp, matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4566 if (aConv->conversionVertex().isValid()) {
4592 float signX = aConv->refittedPairMomentum().x() / fabs(aConv->refittedPairMomentum().x());
4593 float signY = aConv->refittedPairMomentum().y() / fabs(aConv->refittedPairMomentum().y());
4594 float signZ = aConv->refittedPairMomentum().z() / fabs(aConv->refittedPairMomentum().z());
4603 float thetaConv = aConv->refittedPairMomentum().Theta();
4604 float thetaSC = matchingPho->superCluster()->position().theta();
4606 sqrt(matchingPho->superCluster()->position().x() * matchingPho->superCluster()->position().x() +
4607 matchingPho->superCluster()->position().y() * matchingPho->superCluster()->position().y());
4608 float zSC = matchingPho->superCluster()->position().z();
4609 float zPV =
sqrt(rSC * rSC + zSC * zSC) *
sin(thetaConv - thetaSC) /
sin(thetaConv);
4614 if (phoIsInBarrel) {
4617 }
else if (phoIsInEndcap) {
4620 }
else if (phoIsInEndcapP) {
4623 }
else if (phoIsInEndcapM) {
4634 float dPhiTracksAtEcal = -99;
4635 float dEtaTracksAtEcal = -99;
4636 if (!aConv->bcMatchingWithTracks().empty() && aConv->bcMatchingWithTracks()[0].
isNonnull() &&
4637 aConv->bcMatchingWithTracks()[1].
isNonnull()) {
4639 float recoPhi1 = aConv->ecalImpactPosition()[0].phi();
4640 float recoPhi2 = aConv->ecalImpactPosition()[1].phi();
4641 float recoEta1 = aConv->ecalImpactPosition()[0].eta();
4642 float recoEta2 = aConv->ecalImpactPosition()[1].eta();
4649 dPhiTracksAtEcal = recoPhi1 - recoPhi2;
4651 dEtaTracksAtEcal = recoEta1 - recoEta2;
4663 if (phoIsInBarrel) {
4667 if (phoIsInEndcap) {
4673 for (
unsigned int i = 0;
i <
tracks.size();
i++) {
4675 itAss = myAss.find(tfrb.
get());
4676 if (itAss == myAss.end())
4679 float trkProvenance = 3;
4680 if (
tracks[0]->
algoName() ==
"outInEcalSeededConv" &&
tracks[1]->algoName() ==
"outInEcalSeededConv")
4682 if (
tracks[0]->
algoName() ==
"inOutEcalSeededConv" &&
tracks[1]->algoName() ==
"inOutEcalSeededConv")
4686 (
tracks[1]->algoName() ==
"outInEcalSeededConv" &&
tracks[0]->algoName() ==
"inOutEcalSeededConv"))
4704 float simPt =
sqrt(((*itAss).second)->momentum().perp2());
4706 float refPt = -9999.;
4707 float px = 0,
py = 0;
4709 if (aConv->conversionVertex().isValid()) {
4710 reco::Track refTrack = aConv->conversionVertex().refittedTracks()[
i];
4715 float ptres = refPt - simPt;
4717 float pterror = aConv->conversionVertex().refittedTracks()[
i].ptError();
4720 if (trkProvenance == 3)
4731 if (!aConv->bcMatchingWithTracks().empty() && aConv->bcMatchingWithTracks()[
i].
isNonnull())
4733 sqrt(aConv->tracks()[
i]->outerMomentum().Mag2()));
4735 if (phoIsInBarrel) {
4741 if (!aConv->bcMatchingWithTracks().empty() && aConv->bcMatchingWithTracks()[
i].
isNonnull())
4743 sqrt(aConv->tracks()[
i]->outerMomentum().Mag2()));
4745 if (phoIsInEndcap) {
4751 if (!aConv->bcMatchingWithTracks().empty() && aConv->bcMatchingWithTracks()[
i].
isNonnull())
4753 sqrt(aConv->tracks()[
i]->outerMomentum().Mag2()));
4766 if (!atLeastOneRecoTwoTrackConversion) {
4767 for (
unsigned int iConv = 0; iConv < conversionsOneLeg.
size(); iConv++) {
4769 const std::vector<edm::RefToBase<reco::Track> >
tracks = aConv->tracks();
4775 std::map<const reco::Track*, TrackingParticleRef> myAss;
4776 for (
unsigned int i = 0;
i <
tracks.size();
i++) {
4781 float eoverp = photonE /
tracks[0]->p();
4783 if (phoIsInBarrel) {
4789 if (fabs(
mcEta_) <= 1.) {
4799 std::vector<std::pair<RefToBase<reco::Track>,
double> > trackV;
4803 trackV = (std::vector<std::pair<RefToBase<reco::Track>,
double> >)
q[
theConvTP_[0]];
4805 trackV = (std::vector<std::pair<RefToBase<reco::Track>,
double> >)
q[
theConvTP_[1]];
4844 for (
unsigned int iConv = 0; iConv <
conversions.size(); iConv++) {
4846 double like = aConv->MVAout();
4850 const std::vector<edm::RefToBase<reco::Track> >
tracks = aConv->tracks();
4860 bool phoIsInBarrel =
false;
4861 bool phoIsInEndcap =
false;
4862 if (fabs(aPho.
superCluster()->position().eta()) < 1.479) {
4863 phoIsInBarrel =
true;
4865 phoIsInEndcap =
true;
4877 if (aConv->conversionVertex().isValid())
4883 for (std::vector<PhotonMCTruth>::const_iterator mcPho = mcPhotons.begin(); mcPho != mcPhotons.end(); mcPho++) {
4885 float mcPhi = (*mcPho).fourMomentum().phi();
4888 mcEta_ = (*mcPho).fourMomentum().pseudoRapidity();
4896 if (fabs(
mcEta_) > END_HI)
4903 if ((*mcPho).isAConversion() != 1)
4906 (fabs(
mcEta_) > BARL && fabs(
mcEta_) <= END_HI && fabs((*mcPho).vertex().z()) < 210)))
4912 if (fabs(
tp->vx() - (*mcPho).vertex().x()) < 0.0001 && fabs(
tp->vy() - (*mcPho).vertex().y()) < 0.0001 &&
4913 fabs(
tp->vz() - (*mcPho).vertex().z()) < 0.0001) {
4923 std::vector<std::pair<RefToBase<reco::Track>,
double> > trackV1, trackV2;
4925 auto itP1 =
p1.find(tk1);
4926 auto itP2 =
p2.find(tk2);
4927 bool good = (itP1 !=
p1.end()) and (not itP1->val.empty()) and (itP2 !=
p2.end()) and (not itP2->val.empty());
4929 itP1 =
p1.find(tk2);
4930 itP2 =
p2.find(tk1);
4931 good = (itP1 !=
p1.end()) and (not itP1->val.empty()) and (itP2 !=
p2.end()) and (not itP2->val.empty());
4934 std::vector<std::pair<TrackingParticleRef, double> >
const& tp1 = itP1->val;
4935 std::vector<std::pair<TrackingParticleRef, double> >
const& tp2 = itP2->val;
4940 if (
abs(tpr1->pdgId()) == 11 &&
abs(tpr2->pdgId()) == 11) {
4941 if ((tpr1->parentVertex()->sourceTracks_end() - tpr1->parentVertex()->sourceTracks_begin() == 1) &&
4942 (tpr2->parentVertex()->sourceTracks_end() - tpr2->parentVertex()->sourceTracks_begin() == 1)) {
4943 if (tpr1->parentVertex().
key() == tpr2->parentVertex().
key() &&
4944 ((*tpr1->parentVertex()->sourceTracks_begin())->
pdgId() == 22)) {
4957 if (aConv->conversionVertex().isValid())
4964 if (aConv->conversionVertex().isValid()) {
4965 float chi2Prob =
ChiSquaredProbability(aConv->conversionVertex().chi2(), aConv->conversionVertex().ndof());
4966 double convR =
sqrt(aConv->conversionVertex().position().perp2());
4967 double scalar = aConv->conversionVertex().position().x() * aConv->pairMomentum().x() +
4968 aConv->conversionVertex().position().y() * aConv->pairMomentum().y();
4971 convR = -
sqrt(aConv->conversionVertex().position().perp2());
4973 sqrt(aConv->conversionVertex().position().perp2()));
4975 if (!aConv->caloCluster().empty()) {
4978 sqrt(aConv->conversionVertex().position().perp2()));
4979 if (fabs(aConv->caloCluster()[0]->eta()) <= 1.) {
4980 h_convVtxYvsX_->
Fill(aConv->conversionVertex().position().y(), aConv->conversionVertex().position().x());
4985 aConv->conversionVertex().position().x());
4987 aConv->conversionVertex().position().x());
4992 if (fabs(aConv->caloCluster()[0]->eta()) > 1.)
4998 if (phoIsInBarrel) {
5002 if (phoIsInEndcap) {
5014 for (reco::GenJetCollection::const_iterator genJetIter =
genJetCollection.begin();
5019 if (fabs(genJetIter->eta()) > 2.5)
5022 float mcJetPhi = genJetIter->phi();
5025 float mcJetPt = genJetIter->pt();
5031 std::vector<reco::Photon> thePhotons;
5038 float phiPho = aPho.
phi();
5039 float etaPho = aPho.
eta();
5050 matchingPho = *iPho;
5063 bool phoIsInBarrel =
false;
5064 bool phoIsInEndcap =
false;
5065 if (fabs(matchingPho.
superCluster()->position().eta()) < 1.479) {
5066 phoIsInBarrel =
true;
5068 phoIsInEndcap =
true;
5071 if (phoIsInBarrel) {
5074 if (!ecalRecHitHandle.
isValid()) {
5077 edm::LogError(
"PhotonProducer") <<
"Error! Can't get the product " <<
l.module;
5081 }
else if (phoIsInEndcap) {
5084 if (!ecalRecHitHandle.
isValid()) {
5087 edm::LogError(
"PhotonProducer") <<
"Error! Can't get the product " <<
l.module;
5093 float photonE = matchingPho.
energy();
5094 float photonEt = matchingPho.
et();
5095 float r9 = matchingPho.
r9();
5096 float r1 = matchingPho.
r1x5();
5097 float r2 = matchingPho.
r2x5();
5104 std::vector<std::pair<DetId, float> >::const_iterator rhIt;
5106 bool atLeastOneDeadChannel =
false;
5110 for (rhIt = (*bcIt)->hitsAndFractions().begin(); rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) {
5113 if (rhIt->first == (*it).id()) {
5114 if ((*it).recoFlag() == 9) {
5115 atLeastOneDeadChannel =
true;
5123 if (atLeastOneDeadChannel) {
5192 if (phoIsInBarrel) {
5216 }
else if (phoIsInEndcap) {
5245 for (
unsigned int iConv = 0; iConv <
conversions.size(); iConv++) {
5248 const std::vector<edm::RefToBase<reco::Track> >
tracks = aConv->tracks();
5249 double like = aConv->MVAout();
5254 if (!aConv->caloCluster().empty()) {
5259 float eoverp = aConv->EoverP();
5263 float dPhiTracksAtVtx = aConv->dPhiTracksAtVtx();
5266 if (phoIsInBarrel) {
5272 }
else if (phoIsInEndcap) {
5280 if (aConv->conversionVertex().isValid()) {
5281 double convR =
sqrt(aConv->conversionVertex().position().perp2());
5282 double scalar = aConv->conversionVertex().position().x() * aConv->pairMomentum().x() +
5283 aConv->conversionVertex().position().y() * aConv->pairMomentum().y();
5285 convR = -
sqrt(aConv->conversionVertex().position().perp2());
5289 sqrt(aConv->conversionVertex().position().perp2()));
5290 if (!aConv->caloCluster().empty() && fabs(aConv->caloCluster()[0]->eta()) <= 1.) {
5292 aConv->conversionVertex().position().x());
5306 if (!(mcIter->pdgId() == 22))
5308 if (mcIter->mother() !=
nullptr and !(mcIter->mother()->pdgId() == 25))
5310 if (fabs(mcIter->eta()) > 2.5)
5313 float mcPhi = mcIter->phi();
5314 float mcEta = mcIter->eta();
5316 float mcEnergy = mcIter->energy();
5318 double dR = 9999999.;
5319 float minDr = 10000.;
5323 for (
unsigned int ipho = 0; ipho < photonHandle->size(); ipho++) {
5326 double dphi = pho->phi() - mcPhi;
5328 dphi = dphi < 0 ? (CLHEP::twopi) + dphi : dphi - CLHEP::twopi;
5330 double deta = pho->superCluster()->position().eta() - mcEta;
5333 if (
dR < 0.1 &&
dR < minDr) {
5346 bool phoIsInBarrel =
false;
5347 bool phoIsInEndcap =
false;
5349 float phoEta = matchingPho->
superCluster()->position().eta();
5350 if (fabs(phoEta) < 1.479) {
5351 phoIsInBarrel =
true;
5353 phoIsInEndcap =
true;
5356 float photonE = matchingPho->
energy();
5358 float photonEt = matchingPho->
energy() / cosh(matchingPho->
eta());
5361 float r9 = matchingPho->
r9();
5362 float full5x5_r9 = matchingPho->
full5x5_r9();
5363 float r1 = matchingPho->
r1x5();
5364 float r2 = matchingPho->
r2x5();
5382 if ((photonEt > 14 && newhOverE < 0.15) || (photonEt > 10 && photonEt < 14 && chargedHadIso < 10)) {
5414 if (phoIsInBarrel) {
5438 if (phoIsInEndcap) {
Public Types inherited from DQMOneEDAnalyzer<>
1.8.16