3173 using namespace edm;
3177 const float BARL = 1.4442;
3179 const float END_HI = 2.5;
3184 e.getByLabel(
"trackAssociatorByHitsForPhotonValidation", theHitsAssociator);
3188 LogInfo(
"PhotonValidator") <<
"PhotonValidator Analyzing event number: " <<
e.id() <<
" Global Counter " <<
nEvt_
3206 if (!photonHandle.
isValid()) {
3207 edm::LogError(
"PhotonProducer") <<
"Error! Can't get the Photon collection " << std::endl;
3214 if (!pfCandidateHandle.
isValid()) {
3215 edm::LogError(
"PhotonValidator") <<
"Error! Can't get the product pfCandidates " << std::endl;
3220 if (
fName_ ==
"pfPhotonValidator") {
3222 if (!phoToParticleBasedIsoMapHandle.
isValid()) {
3223 edm::LogInfo(
"PhotonValidator") <<
"Error! Can't get the product: valueMap photons to particle based iso "
3226 phoToParticleBasedIsoMap = *(phoToParticleBasedIsoMapHandle.
product());
3244 if ((*itHits)->isValid()) {
3259 if ((*itHits)->isValid()) {
3270 std::vector<SimTrack> theSimTracks;
3271 std::vector<SimVertex> theSimVertices;
3283 theSimTracks.insert(theSimTracks.end(), SimTk->begin(), SimTk->end());
3284 theSimVertices.insert(theSimVertices.end(), SimVtx->begin(), SimVtx->end());
3311 std::vector<reco::PhotonCollection::const_iterator> StoRMatchedConvertedPhotons;
3321 OISimToReco =
trackAssociator->associateSimToReco(outInTrkHandle, ElectronTPHandle);
3322 IOSimToReco =
trackAssociator->associateSimToReco(inOutTrkHandle, ElectronTPHandle);
3324 OIRecoToSim =
trackAssociator->associateRecoToSim(outInTrkHandle, ElectronTPHandle);
3325 IORecoToSim =
trackAssociator->associateRecoToSim(inOutTrkHandle, ElectronTPHandle);
3329 vector<reco::SimToRecoCollection*> StoRCollPtrs;
3330 StoRCollPtrs.push_back(&OISimToReco);
3331 StoRCollPtrs.push_back(&IOSimToReco);
3332 vector<reco::RecoToSimCollection*> RtoSCollPtrs;
3333 RtoSCollPtrs.push_back(&OIRecoToSim);
3334 RtoSCollPtrs.push_back(&IORecoToSim);
3336 for (
int i = 0;
i < 2;
i++)
3338 for (
int i = 0;
i < 2;
i++)
3341 std::vector<reco::PhotonRef> myPhotons;
3343 for (
unsigned int iPho = 0; iPho < photonHandle->size(); iPho++) {
3346 if (fabs(phoRef->eta()) > 2.5)
3348 myPhotons.push_back(phoRef);
3351 std::sort(myPhotons.begin(), myPhotons.end(), sortPhotons());
3353 if (myPhotons.size() >= 2) {
3354 if (myPhotons[0]->
et() > 40 && myPhotons[1]->et() > 25) {
3360 float gamgamMass2 = p12.Dot(p12);
3361 float gamgamMass2_regr1 = p12_regr1.Dot(p12_regr1);
3362 float gamgamMass2_regr2 = p12_regr2.Dot(p12_regr2);
3365 if (gamgamMass2 > 0) {
3368 if (myPhotons[0]->isEB() && myPhotons[1]->isEB())
3370 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3371 (myPhotons[0]->isEB() && myPhotons[1]->isEE()))
3374 if (myPhotons[0]->
r9() > 0.94 && myPhotons[1]->
r9() > 0.94) {
3376 if (myPhotons[0]->isEB() && myPhotons[1]->isEB())
3378 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3379 (myPhotons[0]->isEB() && myPhotons[1]->isEE()))
3389 if (chi2Prob1 > 0.0005 && chi2Prob2 > 0.0005) {
3391 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3394 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3395 (myPhotons[0]->isEB() && myPhotons[1]->isEE())) {
3401 myPhotons[1]->
r9() > 0.93) {
3405 if (chi2Prob1 > 0.0005) {
3407 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3410 if (myPhotons[0]->isEE() || myPhotons[1]->isEE()) {
3416 myPhotons[0]->
r9() > 0.93) {
3420 if (chi2Prob1 > 0.0005) {
3422 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3425 if (myPhotons[0]->isEE() || myPhotons[1]->isEE()) {
3434 if (gamgamMass2_regr1 > 0) {
3437 if (myPhotons[0]->isEB() && myPhotons[1]->isEB())
3439 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3440 (myPhotons[0]->isEB() && myPhotons[1]->isEE()))
3443 if (myPhotons[0]->
r9() > 0.94 && myPhotons[1]->
r9() > 0.94) {
3445 if (myPhotons[0]->isEB() && myPhotons[1]->isEB())
3447 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3448 (myPhotons[0]->isEB() && myPhotons[1]->isEE()))
3459 if (chi2Prob1 > 0.0005 && chi2Prob2 > 0.0005) {
3461 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3464 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3465 (myPhotons[0]->isEB() && myPhotons[1]->isEE())) {
3471 myPhotons[1]->
r9() > 0.93) {
3475 if (chi2Prob1 > 0.0005) {
3477 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3480 if (myPhotons[0]->isEE() || myPhotons[1]->isEE()) {
3486 myPhotons[0]->
r9() > 0.93) {
3490 if (chi2Prob1 > 0.0005) {
3492 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3495 if (myPhotons[0]->isEE() || myPhotons[1]->isEE()) {
3504 if (gamgamMass2_regr2 > 0) {
3507 if (myPhotons[0]->isEB() && myPhotons[1]->isEB())
3509 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3510 (myPhotons[0]->isEB() && myPhotons[1]->isEE()))
3513 if (myPhotons[0]->
r9() > 0.94 && myPhotons[1]->
r9() > 0.94) {
3515 if (myPhotons[0]->isEB() && myPhotons[1]->isEB())
3517 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3518 (myPhotons[0]->isEB() && myPhotons[1]->isEE()))
3529 if (chi2Prob1 > 0.0005 && chi2Prob2 > 0.0005) {
3531 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3534 if ((myPhotons[0]->isEE() && myPhotons[1]->isEE()) || (myPhotons[0]->isEE() && myPhotons[1]->isEB()) ||
3535 (myPhotons[0]->isEB() && myPhotons[1]->isEE())) {
3541 myPhotons[1]->
r9() > 0.93) {
3545 if (chi2Prob1 > 0.0005) {
3547 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3550 if (myPhotons[0]->isEE() || myPhotons[1]->isEE()) {
3556 myPhotons[0]->
r9() > 0.93) {
3560 if (chi2Prob1 > 0.0005) {
3562 if (myPhotons[0]->isEB() && myPhotons[1]->isEB()) {
3565 if (myPhotons[0]->isEE() || myPhotons[1]->isEE()) {
3576 for (std::vector<PhotonMCTruth>::const_iterator mcPho = mcPhotons.begin(); mcPho != mcPhotons.end(); mcPho++) {
3580 for (HepMC::GenEvent::particle_const_iterator mcIter = myGenEvent->particles_begin();
3581 mcIter != myGenEvent->particles_end();
3583 if ((*mcIter)->pdg_id() != 22)
3585 bool isTheSame =
false;
3587 if ((*mcIter)->production_vertex()) {
3588 if ((*mcIter)->production_vertex()->particles_begin(
HepMC::parents) !=
3590 mother = *((*mcIter)->production_vertex()->particles_begin(
HepMC::parents));
3593 float mcPhi = (*mcPho).fourMomentum().phi();
3595 mcEta_ = (*mcPho).fourMomentum().pseudoRapidity();
3597 mcConvR_ = (*mcPho).vertex().perp();
3604 if (fabs(
mcEta_) > END_HI)
3607 if (mother ==
nullptr || (mother !=
nullptr && mother->pdg_id() == 22) ||
3608 (mother !=
nullptr && mother->pdg_id() == 25) || (mother !=
nullptr && mother->pdg_id() == 35)) {
3609 double dPt = fabs((*mcIter)->momentum().perp() - (*mcPho).fourMomentum().et());
3610 float phiMother = (*mcIter)->momentum().phi();
3612 double dEta = fabs((*mcIter)->momentum().eta() - (*mcPho).fourMomentum().pseudoRapidity());
3614 if (
dEta <= 0.0001 &&
dPhi <= 0.0001 && dPt <= 0.0001)
3632 bool goodSimConversion =
false;
3633 bool visibleConversion =
false;
3634 bool visibleConversionsWithTwoSimTracks =
false;
3635 if ((*mcPho).isAConversion() == 1) {
3649 (fabs(
mcEta_) > BARL && fabs(
mcEta_) <= END_HI && fabs((*mcPho).vertex().z()) < 210))
3650 visibleConversion =
true;
3655 if (fabs(
tp->vx() - (*mcPho).vertex().x()) < 0.001 && fabs(
tp->vy() - (*mcPho).vertex().y()) < 0.001 &&
3656 fabs(
tp->vz() - (*mcPho).vertex().z()) < 0.001) {
3661 visibleConversionsWithTwoSimTracks =
true;
3662 goodSimConversion =
false;
3664 if (visibleConversion && visibleConversionsWithTwoSimTracks)
3665 goodSimConversion =
true;
3666 if (goodSimConversion) {
3687 float minDelta = 10000.;
3688 std::vector<reco::PhotonRef> thePhotons;
3693 for (
unsigned int iPho = 0; iPho < photonHandle->size(); iPho++) {
3695 thePhotons.push_back(aPho);
3696 float phiPho = aPho->phi();
3697 float etaPho = aPho->eta();
3730 bool phoIsInBarrel =
false;
3731 bool phoIsInEndcap =
false;
3732 bool phoIsInEndcapP =
false;
3733 bool phoIsInEndcapM =
false;
3737 if (fabs(matchingPho->superCluster()->position().eta()) < 1.479) {
3738 phoIsInBarrel =
true;
3740 phoIsInEndcap =
true;
3741 if (matchingPho->superCluster()->position().eta() > 0)
3742 phoIsInEndcapP =
true;
3743 if (matchingPho->superCluster()->position().eta() < 0)
3744 phoIsInEndcapM =
true;
3748 if (phoIsInBarrel) {
3751 if (!ecalRecHitHandle.
isValid()) {
3754 edm::LogError(
"PhotonProducer") <<
"Error! Can't get the product " <<
l.module;
3758 }
else if (phoIsInEndcap) {
3761 if (!ecalRecHitHandle.
isValid()) {
3764 edm::LogError(
"PhotonProducer") <<
"Error! Can't get the product " <<
l.module;
3771 float photonE = matchingPho->energy();
3772 float sigmaEoE = matchingPho->getCorrectedEnergyError(matchingPho->getCandidateP4type()) / matchingPho->energy();
3774 float photonEt = matchingPho->pt();
3777 float r9 = matchingPho->r9();
3779 float r1 = matchingPho->r1x5();
3780 float r2 = matchingPho->r2x5();
3781 float sigmaIetaIeta = matchingPho->sigmaIetaIeta();
3783 float hOverE = matchingPho->hadronicOverEm();
3784 float newhOverE = matchingPho->hadTowOverEm();
3785 float ecalIso = matchingPho->ecalRecHitSumEtConeDR04();
3786 float hcalIso = matchingPho->hcalTowerSumEtConeDR04();
3787 float newhcalIso = matchingPho->hcalTowerSumEtBcConeDR04();
3788 float trkIso = matchingPho->trkSumPtSolidConeDR04();
3789 float nIsoTrk = matchingPho->nTrkSolidConeDR04();
3791 float chargedHadIso = matchingPho->chargedHadronIso();
3792 float neutralHadIso = matchingPho->neutralHadronIso();
3793 float photonIso = matchingPho->photonIso();
3794 float etOutsideMustache = matchingPho->etOutsideMustache();
3795 int nClusterOutsideMustache = matchingPho->nClusterOutsideMustache();
3796 float pfMVA = matchingPho->pfMVA();
3798 std::vector<std::pair<DetId, float> >::const_iterator rhIt;
3799 bool atLeastOneDeadChannel =
false;
3801 bcIt != matchingPho->superCluster()->clustersEnd();
3803 for (rhIt = (*bcIt)->hitsAndFractions().begin(); rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) {
3806 if (rhIt->first == (*it).id()) {
3807 if ((*it).recoFlag() == 9) {
3808 atLeastOneDeadChannel =
true;
3816 if (atLeastOneDeadChannel) {
3834 h_scEt_[
type][0]->
Fill(matchingPho->superCluster()->energy() / cosh(matchingPho->superCluster()->eta()));
3836 h_psE_->
Fill(matchingPho->superCluster()->preshowerEnergy());
3909 h_phoDEta_[0]->
Fill(matchingPho->eta() - (*mcPho).fourMomentum().eta());
3913 h_nConv_[0][0]->
Fill(
float(matchingPho->conversions().size()));
3914 h_nConv_[1][0]->
Fill(
float(matchingPho->conversionsOneLeg().size()));
3924 p_eResVsEt_[0][0]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
3929 h2_eResVsEt_[0][0]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
3933 h2_sceResVsR9_[0]->
Fill(
r9, matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
3937 p_sceResVsR9_[0]->
Fill(
r9, matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
3939 if ((*mcPho).isAConversion() == 0) {
3970 if (photonE / (*mcPho).fourMomentum().e() < 0.3 && photonE / (*mcPho).fourMomentum().e() > 0.1) {
3973 if ((
r9 > 0.94 && phoIsInBarrel) || (
r9 > 0.95 && phoIsInEndcap)) {
3979 h2_eResVsEt_[0][1]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
3980 p_eResVsEt_[0][1]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
3985 }
else if ((
r9 <= 0.94 && phoIsInBarrel) || (
r9 <= 0.95 && phoIsInEndcap)) {
3990 p_eResVsEt_[0][2]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
3996 h2_eResVsEt_[0][2]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4001 if (phoIsInBarrel) {
4003 h_scEt_[
type][1]->
Fill(matchingPho->superCluster()->energy() / cosh(matchingPho->superCluster()->eta()));
4027 h_nConv_[1][1]->
Fill(
float(matchingPho->conversionsOneLeg().size()));
4033 p_sceResVsR9_[1]->
Fill(
r9, matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4036 h2_sceResVsR9_[1]->
Fill(
r9, matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4039 h2_eResVsEt_[1][0]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4041 p_eResVsEt_[1][0]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4042 p_eResVsNVtx_[1][0]->
Fill(
float(vtxH->size()), photonE / (*mcPho).fourMomentum().e());
4052 h2_eResVsEt_[1][1]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4053 p_eResVsEt_[1][1]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4054 p_eResVsNVtx_[1][1]->
Fill(
float(vtxH->size()), photonE / (*mcPho).fourMomentum().e());
4063 p_eResVsEt_[1][2]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4064 p_eResVsNVtx_[1][2]->
Fill(
float(vtxH->size()), photonE / (*mcPho).fourMomentum().e());
4068 h2_eResVsEt_[1][2]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4073 if (phoIsInEndcap) {
4075 h_scEt_[
type][2]->
Fill(matchingPho->superCluster()->energy() / cosh(matchingPho->superCluster()->eta()));
4099 h_nConv_[1][2]->
Fill(
float(matchingPho->conversionsOneLeg().size()));
4105 p_sceResVsR9_[2]->
Fill(
r9, matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4108 h2_sceResVsR9_[2]->
Fill(
r9, matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4111 h2_eResVsEt_[2][0]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4114 p_eResVsEt_[2][0]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4115 p_eResVsNVtx_[2][0]->
Fill(
float(vtxH->size()), photonE / (*mcPho).fourMomentum().e());
4125 h2_eResVsEt_[2][1]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4126 p_eResVsEt_[2][1]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4127 p_eResVsNVtx_[2][1]->
Fill(
float(vtxH->size()), photonE / (*mcPho).fourMomentum().e());
4136 p_eResVsEt_[2][2]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4137 p_eResVsNVtx_[2][2]->
Fill(
float(vtxH->size()), photonE / (*mcPho).fourMomentum().e());
4142 h2_eResVsEt_[2][2]->
Fill((*mcPho).fourMomentum().et(), photonE / (*mcPho).fourMomentum().e());
4149 if (
fName_ ==
"pfPhotonValidator") {
4150 float SumPtIsoValCh = 0.;
4151 float SumPtIsoValNh = 0.;
4152 float SumPtIsoValPh = 0.;
4154 float SumPtIsoValCleanCh = 0.;
4155 float SumPtIsoValCleanNh = 0.;
4156 float SumPtIsoValCleanPh = 0.;
4158 for (
unsigned int lCand = 0; lCand < pfCandidateHandle->size(); lCand++) {
4160 float dR =
deltaR(matchingPho->eta(), matchingPho->phi(), pfCandRef->eta(), pfCandRef->phi());
4170 SumPtIsoValCh += pfCandRef->pt();
4178 SumPtIsoValNh += pfCandRef->pt();
4186 SumPtIsoValPh += pfCandRef->pt();
4195 for (std::vector<reco::PFCandidateRef>::const_iterator
i = phoToParticleBasedIsoMap[matchingPho].
begin();
4196 i != phoToParticleBasedIsoMap[matchingPho].
end();
4198 if ((*
i) == pfCandRef) {
4206 SumPtIsoValCleanCh += pfCandRef->pt();
4214 SumPtIsoValCleanNh += pfCandRef->pt();
4222 SumPtIsoValCleanPh += pfCandRef->pt();
4239 if (phoIsInBarrel) {
4257 if (!(visibleConversion && visibleConversionsWithTwoSimTracks))
4268 if (fabs(
mcEta_) <= 1.) {
4278 bool atLeastOneRecoTwoTrackConversion =
false;
4279 for (
unsigned int iConv = 0; iConv <
conversions.size(); iConv++) {
4281 double like = aConv->MVAout();
4290 const std::vector<edm::RefToBase<reco::Track> >
tracks = aConv->tracks();
4293 atLeastOneRecoTwoTrackConversion =
true;
4297 if (
tracks.size() == 2) {
4312 std::map<const reco::Track*, TrackingParticleRef> myAss;
4313 std::map<const reco::Track*, TrackingParticleRef>::const_iterator itAss;
4314 std::map<reco::TrackRef, TrackingParticleRef>::const_iterator itAssMin;
4315 std::map<reco::TrackRef, TrackingParticleRef>::const_iterator itAssMax;
4319 for (
unsigned int i = 0;
i <
tracks.size();
i++) {
4336 std::vector<std::pair<RefToBase<reco::Track>,
double> > trackV;
4340 trackV = (std::vector<std::pair<RefToBase<reco::Track>,
double> >)
q[
theConvTP_[0]];
4342 trackV = (std::vector<std::pair<RefToBase<reco::Track>,
double> >)
q[
theConvTP_[1]];
4356 float refP = -99999.;
4357 float refPt = -99999.;
4358 if (aConv->conversionVertex().isValid()) {
4359 refP =
sqrt(aConv->refittedPairMomentum().Mag2());
4360 refPt =
sqrt(aConv->refittedPairMomentum().perp2());
4362 float invM = aConv->pairInvariantMass();
4371 if (
tracks.size() == 2) {
4378 if (!aConv->caloCluster().empty())
4381 float trkProvenance = 3;
4382 if (
tracks[0]->
algoName() ==
"outInEcalSeededConv" &&
tracks[1]->algoName() ==
"outInEcalSeededConv")
4384 if (
tracks[0]->
algoName() ==
"inOutEcalSeededConv" &&
tracks[1]->algoName() ==
"inOutEcalSeededConv")
4387 (
tracks[1]->algoName() ==
"outInEcalSeededConv" &&
tracks[0]->algoName() ==
"inOutEcalSeededConv"))
4389 if (trkProvenance == 3) {
4406 if (!aConv->caloCluster().empty())
4416 if (aConv->conversionVertex().isValid()) {
4426 if (chi2Prob > 0.0005) {
4434 if (chi2Prob > 0.0005) {
4435 if (!aConv->caloCluster().empty()) {
4438 h_convERes_[0][0]->
Fill(aConv->caloCluster()[0]->energy() / (*mcPho).fourMomentum().e());
4448 if (phoIsInBarrel) {
4449 if (!aConv->caloCluster().empty())
4450 h_convERes_[0][1]->
Fill(aConv->caloCluster()[0]->energy() / (*mcPho).fourMomentum().e());
4455 if (phoIsInEndcap) {
4456 if (!aConv->caloCluster().empty())
4457 h_convERes_[0][2]->
Fill(aConv->caloCluster()[0]->energy() / (*mcPho).fourMomentum().e());
4471 float eoverp = -99999.;
4473 if (aConv->conversionVertex().isValid()) {
4474 eoverp = photonE /
sqrt(aConv->refittedPairMomentum().Mag2());
4481 matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4498 matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4501 matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4509 float dPhiTracksAtVtx = aConv->dPhiTracksAtVtx();
4526 if (phoIsInBarrel) {
4528 if (aConv->conversionVertex().isValid()) {
4540 eoverp, matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4548 if (phoIsInEndcap) {
4550 if (aConv->conversionVertex().isValid()) {
4556 eoverp, matchingPho->superCluster()->energy() / (*mcPho).fourMomentum().e());
4564 if (aConv->conversionVertex().isValid()) {
4590 float signX = aConv->refittedPairMomentum().x() / fabs(aConv->refittedPairMomentum().x());
4591 float signY = aConv->refittedPairMomentum().y() / fabs(aConv->refittedPairMomentum().y());
4592 float signZ = aConv->refittedPairMomentum().z() / fabs(aConv->refittedPairMomentum().z());
4601 float thetaConv = aConv->refittedPairMomentum().Theta();
4602 float thetaSC = matchingPho->superCluster()->position().theta();
4604 sqrt(matchingPho->superCluster()->position().x() * matchingPho->superCluster()->position().x() +
4605 matchingPho->superCluster()->position().y() * matchingPho->superCluster()->position().y());
4606 float zSC = matchingPho->superCluster()->position().z();
4607 float zPV =
sqrt(rSC * rSC + zSC * zSC) *
sin(thetaConv - thetaSC) /
sin(thetaConv);
4612 if (phoIsInBarrel) {
4615 }
else if (phoIsInEndcap) {
4618 }
else if (phoIsInEndcapP) {
4621 }
else if (phoIsInEndcapM) {
4632 float dPhiTracksAtEcal = -99;
4633 float dEtaTracksAtEcal = -99;
4634 if (!aConv->bcMatchingWithTracks().empty() && aConv->bcMatchingWithTracks()[0].
isNonnull() &&
4635 aConv->bcMatchingWithTracks()[1].
isNonnull()) {
4637 float recoPhi1 = aConv->ecalImpactPosition()[0].phi();
4638 float recoPhi2 = aConv->ecalImpactPosition()[1].phi();
4639 float recoEta1 = aConv->ecalImpactPosition()[0].eta();
4640 float recoEta2 = aConv->ecalImpactPosition()[1].eta();
4647 dPhiTracksAtEcal = recoPhi1 - recoPhi2;
4649 dEtaTracksAtEcal = recoEta1 - recoEta2;
4661 if (phoIsInBarrel) {
4665 if (phoIsInEndcap) {
4671 for (
unsigned int i = 0;
i <
tracks.size();
i++) {
4673 itAss = myAss.find(tfrb.
get());
4674 if (itAss == myAss.end())
4677 float trkProvenance = 3;
4678 if (
tracks[0]->
algoName() ==
"outInEcalSeededConv" &&
tracks[1]->algoName() ==
"outInEcalSeededConv")
4680 if (
tracks[0]->
algoName() ==
"inOutEcalSeededConv" &&
tracks[1]->algoName() ==
"inOutEcalSeededConv")
4684 (
tracks[1]->algoName() ==
"outInEcalSeededConv" &&
tracks[0]->algoName() ==
"inOutEcalSeededConv"))
4702 float simPt =
sqrt(((*itAss).second)->momentum().perp2());
4704 float refPt = -9999.;
4705 float px = 0,
py = 0;
4707 if (aConv->conversionVertex().isValid()) {
4708 reco::Track refTrack = aConv->conversionVertex().refittedTracks()[
i];
4713 float ptres = refPt - simPt;
4715 float pterror = aConv->conversionVertex().refittedTracks()[
i].ptError();
4718 if (trkProvenance == 3)
4729 if (!aConv->bcMatchingWithTracks().empty() && aConv->bcMatchingWithTracks()[
i].
isNonnull())
4731 sqrt(aConv->tracks()[
i]->outerMomentum().Mag2()));
4733 if (phoIsInBarrel) {
4739 if (!aConv->bcMatchingWithTracks().empty() && aConv->bcMatchingWithTracks()[
i].
isNonnull())
4741 sqrt(aConv->tracks()[
i]->outerMomentum().Mag2()));
4743 if (phoIsInEndcap) {
4749 if (!aConv->bcMatchingWithTracks().empty() && aConv->bcMatchingWithTracks()[
i].
isNonnull())
4751 sqrt(aConv->tracks()[
i]->outerMomentum().Mag2()));
4764 if (!atLeastOneRecoTwoTrackConversion) {
4765 for (
unsigned int iConv = 0; iConv < conversionsOneLeg.
size(); iConv++) {
4767 const std::vector<edm::RefToBase<reco::Track> >
tracks = aConv->tracks();
4773 std::map<const reco::Track*, TrackingParticleRef> myAss;
4774 for (
unsigned int i = 0;
i <
tracks.size();
i++) {
4779 float eoverp = photonE /
tracks[0]->p();
4781 if (phoIsInBarrel) {
4787 if (fabs(
mcEta_) <= 1.) {
4797 std::vector<std::pair<RefToBase<reco::Track>,
double> > trackV;
4801 trackV = (std::vector<std::pair<RefToBase<reco::Track>,
double> >)
q[
theConvTP_[0]];
4803 trackV = (std::vector<std::pair<RefToBase<reco::Track>,
double> >)
q[
theConvTP_[1]];
4842 for (
unsigned int iConv = 0; iConv <
conversions.size(); iConv++) {
4844 double like = aConv->MVAout();
4848 const std::vector<edm::RefToBase<reco::Track> >
tracks = aConv->tracks();
4858 bool phoIsInBarrel =
false;
4859 bool phoIsInEndcap =
false;
4860 if (fabs(aPho.
superCluster()->position().eta()) < 1.479) {
4861 phoIsInBarrel =
true;
4863 phoIsInEndcap =
true;
4875 if (aConv->conversionVertex().isValid())
4881 for (std::vector<PhotonMCTruth>::const_iterator mcPho = mcPhotons.begin(); mcPho != mcPhotons.end(); mcPho++) {
4883 float mcPhi = (*mcPho).fourMomentum().phi();
4886 mcEta_ = (*mcPho).fourMomentum().pseudoRapidity();
4894 if (fabs(
mcEta_) > END_HI)
4901 if ((*mcPho).isAConversion() != 1)
4904 (fabs(
mcEta_) > BARL && fabs(
mcEta_) <= END_HI && fabs((*mcPho).vertex().z()) < 210)))
4910 if (fabs(
tp->vx() - (*mcPho).vertex().x()) < 0.0001 && fabs(
tp->vy() - (*mcPho).vertex().y()) < 0.0001 &&
4911 fabs(
tp->vz() - (*mcPho).vertex().z()) < 0.0001) {
4921 std::vector<std::pair<RefToBase<reco::Track>,
double> > trackV1, trackV2;
4923 auto itP1 =
p1.find(tk1);
4924 auto itP2 =
p2.find(tk2);
4925 bool good = (itP1 !=
p1.end()) and (not itP1->val.empty()) and (itP2 !=
p2.end()) and (not itP2->val.empty());
4927 itP1 =
p1.find(tk2);
4928 itP2 =
p2.find(tk1);
4929 good = (itP1 !=
p1.end()) and (not itP1->val.empty()) and (itP2 !=
p2.end()) and (not itP2->val.empty());
4932 std::vector<std::pair<TrackingParticleRef, double> >
const& tp1 = itP1->val;
4933 std::vector<std::pair<TrackingParticleRef, double> >
const& tp2 = itP2->val;
4938 if (
abs(tpr1->pdgId()) == 11 &&
abs(tpr2->pdgId()) == 11) {
4939 if ((tpr1->parentVertex()->sourceTracks_end() - tpr1->parentVertex()->sourceTracks_begin() == 1) &&
4940 (tpr2->parentVertex()->sourceTracks_end() - tpr2->parentVertex()->sourceTracks_begin() == 1)) {
4941 if (tpr1->parentVertex().
key() == tpr2->parentVertex().
key() &&
4942 ((*tpr1->parentVertex()->sourceTracks_begin())->
pdgId() == 22)) {
4955 if (aConv->conversionVertex().isValid())
4962 if (aConv->conversionVertex().isValid()) {
4963 float chi2Prob =
ChiSquaredProbability(aConv->conversionVertex().chi2(), aConv->conversionVertex().ndof());
4964 double convR =
sqrt(aConv->conversionVertex().position().perp2());
4965 double scalar = aConv->conversionVertex().position().x() * aConv->pairMomentum().x() +
4966 aConv->conversionVertex().position().y() * aConv->pairMomentum().y();
4969 convR = -
sqrt(aConv->conversionVertex().position().perp2());
4971 sqrt(aConv->conversionVertex().position().perp2()));
4973 if (!aConv->caloCluster().empty()) {
4976 sqrt(aConv->conversionVertex().position().perp2()));
4977 if (fabs(aConv->caloCluster()[0]->eta()) <= 1.) {
4978 h_convVtxYvsX_->
Fill(aConv->conversionVertex().position().y(), aConv->conversionVertex().position().x());
4983 aConv->conversionVertex().position().x());
4985 aConv->conversionVertex().position().x());
4990 if (fabs(aConv->caloCluster()[0]->eta()) > 1.)
4996 if (phoIsInBarrel) {
5000 if (phoIsInEndcap) {
5012 for (reco::GenJetCollection::const_iterator genJetIter =
genJetCollection.begin();
5017 if (fabs(genJetIter->eta()) > 2.5)
5020 float mcJetPhi = genJetIter->phi();
5023 float mcJetPt = genJetIter->pt();
5029 std::vector<reco::Photon> thePhotons;
5036 float phiPho = aPho.
phi();
5037 float etaPho = aPho.
eta();
5048 matchingPho = *iPho;
5061 bool phoIsInBarrel =
false;
5062 bool phoIsInEndcap =
false;
5063 if (fabs(matchingPho.
superCluster()->position().eta()) < 1.479) {
5064 phoIsInBarrel =
true;
5066 phoIsInEndcap =
true;
5069 if (phoIsInBarrel) {
5072 if (!ecalRecHitHandle.
isValid()) {
5075 edm::LogError(
"PhotonProducer") <<
"Error! Can't get the product " <<
l.module;
5079 }
else if (phoIsInEndcap) {
5082 if (!ecalRecHitHandle.
isValid()) {
5085 edm::LogError(
"PhotonProducer") <<
"Error! Can't get the product " <<
l.module;
5091 float photonE = matchingPho.
energy();
5092 float photonEt = matchingPho.
et();
5093 float r9 = matchingPho.
r9();
5094 float r1 = matchingPho.
r1x5();
5095 float r2 = matchingPho.
r2x5();
5102 std::vector<std::pair<DetId, float> >::const_iterator rhIt;
5104 bool atLeastOneDeadChannel =
false;
5108 for (rhIt = (*bcIt)->hitsAndFractions().begin(); rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) {
5111 if (rhIt->first == (*it).id()) {
5112 if ((*it).recoFlag() == 9) {
5113 atLeastOneDeadChannel =
true;
5121 if (atLeastOneDeadChannel) {
5190 if (phoIsInBarrel) {
5214 }
else if (phoIsInEndcap) {
5243 for (
unsigned int iConv = 0; iConv <
conversions.size(); iConv++) {
5246 const std::vector<edm::RefToBase<reco::Track> >
tracks = aConv->tracks();
5247 double like = aConv->MVAout();
5252 if (!aConv->caloCluster().empty()) {
5257 float eoverp = aConv->EoverP();
5261 float dPhiTracksAtVtx = aConv->dPhiTracksAtVtx();
5264 if (phoIsInBarrel) {
5270 }
else if (phoIsInEndcap) {
5278 if (aConv->conversionVertex().isValid()) {
5279 double convR =
sqrt(aConv->conversionVertex().position().perp2());
5280 double scalar = aConv->conversionVertex().position().x() * aConv->pairMomentum().x() +
5281 aConv->conversionVertex().position().y() * aConv->pairMomentum().y();
5283 convR = -
sqrt(aConv->conversionVertex().position().perp2());
5287 sqrt(aConv->conversionVertex().position().perp2()));
5288 if (!aConv->caloCluster().empty() && fabs(aConv->caloCluster()[0]->eta()) <= 1.) {
5290 aConv->conversionVertex().position().x());
5304 if (!(mcIter->pdgId() == 22))
5306 if (mcIter->mother() !=
nullptr and !(mcIter->mother()->pdgId() == 25))
5308 if (fabs(mcIter->eta()) > 2.5)
5311 float mcPhi = mcIter->phi();
5312 float mcEta = mcIter->eta();
5314 float mcEnergy = mcIter->energy();
5316 double dR = 9999999.;
5317 float minDr = 10000.;
5321 for (
unsigned int ipho = 0; ipho < photonHandle->size(); ipho++) {
5324 double dphi = pho->phi() - mcPhi;
5326 dphi = dphi < 0 ? (CLHEP::twopi) + dphi : dphi - CLHEP::twopi;
5328 double deta = pho->superCluster()->position().eta() - mcEta;
5331 if (
dR < 0.1 &&
dR < minDr) {
5344 bool phoIsInBarrel =
false;
5345 bool phoIsInEndcap =
false;
5347 float phoEta = matchingPho->
superCluster()->position().eta();
5348 if (fabs(phoEta) < 1.479) {
5349 phoIsInBarrel =
true;
5351 phoIsInEndcap =
true;
5354 float photonE = matchingPho->
energy();
5356 float photonEt = matchingPho->
energy() / cosh(matchingPho->
eta());
5359 float r9 = matchingPho->
r9();
5360 float full5x5_r9 = matchingPho->
full5x5_r9();
5361 float r1 = matchingPho->
r1x5();
5362 float r2 = matchingPho->
r2x5();
5380 if ((photonEt > 14 && newhOverE < 0.15) || (photonEt > 10 && photonEt < 14 && chargedHadIso < 10)) {
5412 if (phoIsInBarrel) {
5436 if (phoIsInEndcap) {
Public Types inherited from DQMOneEDAnalyzer<>
1.8.16