32 double CElecPtMin = ps.
getParameter<
double>(
"CElecPtMin");
33 double CEB_siEiE = ps.
getParameter<
double>(
"CEB_sigmaIEtaIEta");
34 double CEB_dPhiIn = ps.
getParameter<
double>(
"CEB_deltaPhiIn");
35 double CEB_dEtaIn = ps.
getParameter<
double>(
"CEB_deltaEtaIn");
36 double CEB_EcalIso = ps.
getParameter<
double>(
"CEB_EcalIso");
37 double CEB_HcalIso = ps.
getParameter<
double>(
"CEB_HcalIso");
38 double CEB_TrckIso = ps.
getParameter<
double>(
"CEB_TrckIso");
39 double CEE_siEiE = ps.
getParameter<
double>(
"CEE_sigmaIEtaIEta");
40 double CEE_dPhiIn = ps.
getParameter<
double>(
"CEE_deltaPhiIn");
41 double CEE_dEtaIn = ps.
getParameter<
double>(
"CEE_deltaEtaIn");
42 double CEE_EcalIso = ps.
getParameter<
double>(
"CEE_EcalIso");
43 double CEE_HcalIso = ps.
getParameter<
double>(
"CEE_HcalIso");
44 double CEE_TrckIso = ps.
getParameter<
double>(
"CEE_TrckIso");
61 for (std::vector<double>::const_iterator it =
CutVector_.begin(); it !=
CutVector_.end(); ++it) {
72 t_ = fs->
make<TTree>(
"ErsatzMEt",
"Data on ErsatzMEt");
76 t_->Branch(
"nTags", &
nTags_,
"nTags/I");
101 t_->Branch(
"CaloMEt", &
CaloMEt_,
"CaloMEt/D");
103 t_->Branch(
"T1MEt", &
T1MEt_,
"T1MEt/D");
105 t_->Branch(
"PfMEt", &
PfMEt_,
"PfMEt/D");
107 t_->Branch(
"TcMEt", &
TcMEt_,
"TcMEt/D");
111 t_->Branch(
"tag_q",
tag_q_,
"tag_q[4]/I");
112 t_->Branch(
"tag_pt",
tag_pt_,
"tag_pt[4]/D");
113 t_->Branch(
"tag_eta",
tag_eta_,
"tag_eta[4]/D");
114 t_->Branch(
"tag_phi",
tag_phi_,
"tag_phi[4]/D");
123 t_->Branch(
"tag_e5x5",
tag_e5x5_,
"tag_e5x5[4]/D");
124 t_->Branch(
"tag_hoe",
tag_hoe_,
"tag_hoe[4]/D");
125 t_->Branch(
"tag_eop",
tag_eop_,
"tag_eop[4]/D");
126 t_->Branch(
"tag_pin",
tag_pin_,
"tag_pin[4]/D");
127 t_->Branch(
"tag_pout",
tag_pout_,
"tag_pout[4]/D");
132 edm::LogDebug_(
"",
"", 103) <<
"Creating ersatz neutrino branches.";
133 t_->Branch(
"probe_q",
probe_q_,
"probe_q[4]/I");
134 t_->Branch(
"probe_pt",
probe_pt_,
"probe_pt[4]/D");
154 t_->Branch(
"Z_m",
Z_m_,
"Z_m[4]/D");
155 t_->Branch(
"Z_pt",
Z_pt_,
"Z_pt[4]/D");
156 t_->Branch(
"Z_eta",
Z_eta_,
"Z_eta[4]/D");
157 t_->Branch(
"Z_y",
Z_y_,
"Z_y[4]/D");
158 t_->Branch(
"Z_phi",
Z_phi_,
"Z_phi[4]/D");
159 t_->Branch(
"Z_rescM",
Z_rescM_,
"Z_rescM[4]/D");
160 t_->Branch(
"Z_rescPt",
Z_rescPt_,
"Z_rescPt[4]/D");
162 t_->Branch(
"Z_rescY",
Z_rescY_,
"Z_rescY[4]/D");
183 t_->Branch(
"McZ_m", &
McZ_m_,
"McZ_m/D");
185 t_->Branch(
"McZ_Pt", &
McZ_pt_,
"McZ_Pt/D");
186 t_->Branch(
"McZ_y", &
McZ_y_,
"McZ_y/D");
187 t_->Branch(
"McZ_Eta", &
McZ_eta_,
"McZ_Eta/D");
188 t_->Branch(
"McZ_Phi", &
McZ_phi_,
"McZ_Phi/D");
354 edm::LogDebug_(
"",
"", 180) <<
"Initialisation of array index " <<
i <<
" completed.";
380 std::vector<math::XYZTLorentzVector> McElecs, McElecsFinalState;
381 std::vector<math::XYZTLorentzVector> McElecsResc;
386 for (reco::GenParticleCollection::const_iterator McP = McCand->begin(); McP != McCand->end(); ++McP) {
389 McElecs.push_back(McP->p4());
393 std::cout <<
"Found electron, ID = " << McP->pdgId() <<
"\t status = " << McP->status() << std::endl;
394 if (McP->status() != 1) {
397 while (McPD->
status() != 1) {
401 for (
int j = 0;
j <
n; ++
j) {
403 std::cout <<
"Daughter " <<
j <<
"\t id = " <<
d->pdgId() << std::endl;
410 std::cout << McPD->
pdgId() <<
" : status = " << McPD->
status() <<
"\tAdding to vector!" << std::endl;
411 McElecsFinalState.push_back(McPD->
p4());
413 McElecsFinalState.push_back(McP->p4());
420 McZ_y_ = Zboson.Rapidity();
426 McElecsResc.resize(2);
428 RescZboson.SetCoordinates(
435 ROOT::Math::Boost CoMBoost(Zboson.BoostToCM());
442 double E_W = RescZboson.E();
443 ROOT::Math::Boost BackToLab(RescZboson.Px() / E_W, RescZboson.Py() / E_W, RescZboson.Pz() / E_W);
445 RescMcElec0 = BackToLab(RescMcElec0);
450 RescMcElec1 = BackToLab(RescMcElec1);
454 McElecsResc[0] = RescMcElec0;
455 McElecsResc[1] = RescMcElec1;
457 edm::LogDebug_(
"",
"", 307) <<
"McElecsResc[0] + McElecsResc[1] = (" << sum.Px() <<
", " << sum.Py() <<
", "
458 << sum.Pz() <<
", " << sum.E() <<
")";
464 for (
unsigned int itrig = 0; itrig < HltRes->
size(); ++itrig) {
466 edm::LogInfo(
"") << itrig <<
" : Name = " << nom <<
"\t Accepted = " << HltRes->
accept(itrig);
469 if (HltRes->
accept(34) == 0)
471 if (HltRes->
accept(34) != 0) {
472 std::vector<reco::GsfElectronRef> UniqueElectrons;
474 edm::LogDebug_(
"",
"ErsatzMEt", 192) <<
"Unique electron size = " << UniqueElectrons.size();
475 std::vector<reco::GsfElectronRef> SelectedElectrons;
477 std::cout <<
"Filter Id = " << fId << std::endl;
479 nTags_ = SelectedElectrons.size();
486 for (std::vector<reco::GsfElectronRef>::const_iterator elec = SelectedElectrons.begin();
487 elec != SelectedElectrons.end();
489 for (
int m = 0;
m < 2; ++
m) {
490 double dRLimit = 99.;
503 std::map<reco::GsfElectronRef, reco::GsfElectronRef>
TagProbePairs;
510 const reco::MET calomet = *(caloMEtCollection->begin());
520 const reco::MET pfmet = *(pfMEtCollection->begin());
525 const reco::MET tcmet = *(tcMEtCollection->begin());
531 for (std::map<reco::GsfElectronRef, reco::GsfElectronRef>::const_iterator it =
TagProbePairs.begin();
545 it->first->isolationVariables04().hcalDepth2TowerSumEt;
577 it->second->isolationVariables04().hcalDepth2TowerSumEt;
600 double dRLimit = 0.2;
601 for (
unsigned int mcEId = 0; mcEId < McElecs.size(); ++mcEId) {
688 std::map<reco::GsfElectronRef, reco::GsfElectronRef> TagProbes;
689 for (std::vector<reco::GsfElectronRef>::const_iterator tagelec =
tags.begin(); tagelec !=
tags.end(); ++tagelec) {
691 std::pair<reco::GsfElectronRef, reco::GsfElectronRef> TagProbePair;
692 int nProbesPerTag = 0;
694 for (reco::GsfElectronCollection::const_iterator probeelec = probeCands->begin(); probeelec != probeCands->end();
697 double probeScEta = probe->superCluster()->eta();
698 if (probe->superCluster() !=
tag->superCluster() && fabs(probeScEta) < 2.5) {
699 if (fabs(probeScEta) < 1.4442 || fabs(probeScEta) > 1.560) {
700 double invmass = ROOT::Math::VectorUtil::InvariantMass(
tag->p4(), probe->p4());
702 TagProbePair = std::make_pair(
tag, probe);
710 if (nProbesPerTag == 1)
711 TagProbes.insert(TagProbePair);
721 elec->TrackPositionAtVtx().X(), elec->TrackPositionAtVtx().Y(), elec->TrackPositionAtVtx().Z());
731 edm::LogDebug_(
"ersatzFabrikV1",
"", 569) <<
"elec = (" << elec->p4().Px() <<
", " << elec->p4().Py() <<
", "
732 << elec->p4().Pz() <<
", " << elec->p4().E() <<
")";
734 edm::LogDebug_(
"ersatzFabrikV1",
"", 569) <<
"Z pt = " << Zboson.Pt() <<
"Z boson mass = " << Zboson.M();
735 edm::LogDebug_(
"ersatzFabrikV1",
"", 570) <<
"Z boson in lab frame = (" << Zboson.Px() <<
", " << Zboson.Py() <<
", "
736 << Zboson.Pz() <<
", " << Zboson.E() <<
")";
738 Zboson.Px(), Zboson.Py(), Zboson.Pz(),
sqrt(Zboson.P2() + (
mW_ *
mW_ * Zboson.M2()) / (
mZ_ *
mZ_)));
739 edm::LogDebug_(
"ersatzFabrikV1",
"", 573) <<
"W boson in lab frame = (" << RescZboson.Px() <<
", " << RescZboson.Py()
740 <<
", " << RescZboson.Pz() <<
", " << RescZboson.E() <<
")";
741 ROOT::Math::Boost BoostToZRestFrame(Zboson.BoostToCM());
742 edm::LogDebug_(
"ersatzFabrikV1",
"", 576) <<
"Electron in lab frame = (" << boost_ele.Px() <<
", " << boost_ele.Py()
743 <<
", " << boost_ele.Pz() <<
", " << boost_ele.E() <<
")";
744 edm::LogDebug_(
"ersatzFabrikV1",
"", 578) <<
"Ersatz Neutrino in lab frame = (" << boost_nu.Px() <<
", "
745 << boost_nu.Py() <<
", " << boost_nu.Pz() <<
", " << boost_nu.E() <<
")";
746 boost_ele = BoostToZRestFrame(boost_ele);
747 boost_nu = BoostToZRestFrame(boost_nu);
748 edm::LogDebug_(
"ersatzFabrikV1",
"", 582) <<
"Electron in Z rest frame = (" << boost_ele.Px() <<
", "
749 << boost_ele.Py() <<
", " << boost_ele.Pz() <<
", " << boost_ele.E() <<
")";
750 edm::LogDebug_(
"ersatzFabrikV1",
"", 584) <<
"Ersatz Neutrino in Z rest frame = (" << boost_nu.Px() <<
", "
751 << boost_nu.Py() <<
", " << boost_nu.Pz() <<
", " << boost_nu.E() <<
")";
755 double E_W = RescZboson.E();
756 ROOT::Math::Boost BackToLab(RescZboson.Px() / E_W, RescZboson.Py() / E_W, RescZboson.Pz() / E_W);
758 boost_ele = BackToLab(boost_ele);
760 boost_nu = BackToLab(boost_nu);
762 edm::LogDebug_(
"ersatzFabrikV1",
"", 597) <<
"Electron back in lab frame = (" << boost_ele.Px() <<
", "
763 << boost_ele.Py() <<
", " << boost_ele.Pz() <<
", " << boost_ele.E() <<
")";
764 edm::LogDebug_(
"ersatzFabrikV1",
"", 599) <<
"Ersatz Neutrino back in lab frame = (" << boost_nu.Px() <<
", "
765 << boost_nu.Py() <<
", " << boost_nu.Pz() <<
", " << boost_nu.E() <<
")";
767 <<
"boost_ele + boost_nu = (" << sum.Px() <<
", " << sum.Py() <<
", " << sum.Pz() <<
", " << sum.E() <<
")";
769 nu.SetXYZT(nu.X(), nu.Y(), 0., nu.T());
770 ele.SetXYZT(ele.X(), ele.Y(), 0., ele.T());
771 boost_ele.SetXYZT(boost_ele.X(), boost_ele.Y(), 0., boost_ele.T());
772 metVec =
met.p4() + nu + ele - boost_ele;
796 sqrt(2. * boost_ele.Pt() * ersatzMEt.pt() * (1 -
cos(
reco::deltaPhi(boost_ele.Phi(), ersatzMEt.phi()))));
819 sqrt(2. * boost_ele.Pt() * ersatzMEt.pt() * (1 -
cos(
reco::deltaPhi(boost_ele.Phi(), ersatzMEt.phi()))));
842 sqrt(2. * boost_ele.Pt() * ersatzMEt.pt() * (1 -
cos(
reco::deltaPhi(boost_ele.Phi(), ersatzMEt.phi()))));
865 sqrt(2. * boost_ele.Pt() * ersatzMEt.pt() * (1 -
cos(
reco::deltaPhi(boost_ele.Phi(), ersatzMEt.phi()))));
874 boost_elec =
tag->p4();
875 edm::LogDebug_(
"ersatzFabrikV1",
"", 858) <<
"boost_elec = (" << boost_elec.Px() <<
", " << boost_elec.Py() <<
", "
876 << boost_elec.Pz() <<
", " << boost_elec.E() <<
")";
877 boost_nu = probe->p4();
878 edm::LogDebug_(
"ersatzFabrikV1",
"", 860) <<
"boost_nu = (" << boost_nu.Px() <<
", " << boost_nu.Py() <<
", "
879 << boost_nu.Pz() <<
", " << boost_nu.E() <<
")";
882 <<
"Zboson = (" << Zboson.Px() <<
", " << Zboson.Py() <<
", " << Zboson.Pz() <<
", " << Zboson.E() <<
")";
884 Zboson.Px(), Zboson.Py(), Zboson.Pz(),
sqrt(Zboson.P2() + (
mW_ *
mW_ * Zboson.M2()) / (
mZ_ *
mZ_)));
885 edm::LogDebug_(
"ersatzFabrikV1",
"", 864) <<
"RescZboson = (" << RescZboson.Px() <<
", " << RescZboson.Py() <<
", "
886 << RescZboson.Pz() <<
", " << RescZboson.E() <<
")";
887 ROOT::Math::Boost BoostToZRestFrame(Zboson.BoostToCM());
888 elec = BoostToZRestFrame(boost_elec);
889 edm::LogDebug_(
"ersatzFabrikV1",
"", 867) <<
"boost_elec (in Z rest frame) = (" << elec.Px() <<
", " << elec.Py()
890 <<
", " << elec.Pz() <<
", " << elec.E() <<
")";
891 nu = BoostToZRestFrame(boost_nu);
893 <<
"boost_nu (in Z rest frame) = (" << nu.Px() <<
", " << nu.Py() <<
", " << nu.Pz() <<
", " << nu.E() <<
")";
896 <<
"elec (in Z rest frame) = (" << elec.Px() <<
", " << elec.Py() <<
", " << elec.Pz() <<
", " << elec.E() <<
")";
899 <<
"nu (in Z rest frame) = (" << nu.Px() <<
", " << nu.Py() <<
", " << nu.Pz() <<
", " << nu.E() <<
")";
900 ROOT::Math::Boost BoostBackToLab(
901 RescZboson.Px() / RescZboson.E(), RescZboson.Py() / RescZboson.E(), RescZboson.Pz() / RescZboson.E());
903 elec = BoostBackToLab(elec);
905 <<
"elec = (" << elec.Px() <<
", " << elec.Py() <<
", " << elec.Pz() <<
", " << elec.E() <<
")";
906 nu = BoostBackToLab(nu);
908 <<
"nu = (" << nu.Px() <<
", " << nu.Py() <<
", " << nu.Pz() <<
", " << nu.E() <<
")";
926 elec.SetXYZT(elec.X(), elec.Y(), 0., elec.T());
927 nu.SetXYZT(nu.X(), nu.Y(), 0., nu.T());
928 boost_elec.SetXYZT(boost_elec.X(), boost_elec.Y(), 0., boost_elec.T());
929 metVec =
met.p4() + nu + elec - boost_elec;