PageRenderTime 177ms CodeModel.GetById 14ms app.highlight 147ms RepoModel.GetById 1ms app.codeStats 1ms

/ElectroWeakAnalysis/ZEE/src/ErsatzMEt.cc

https://github.com/aivanov-cern/cmssw
C++ | 948 lines | 777 code | 68 blank | 103 comment | 60 complexity | 2313a34fb4c823d43c26c4f86f2bc9f0 MD5 | raw file
  1#include "ElectroWeakAnalysis/ZEE/interface/ErsatzMEt.h"
  2#include "FWCore/Common/interface/TriggerNames.h"
  3
  4ErsatzMEt::ErsatzMEt(const edm::ParameterSet& ps)
  5{
  6	MCTruthCollection_  = ps.getParameter<edm::InputTag>("MCTruthCollection");
  7	ElectronCollection_ = ps.getParameter<edm::InputTag>("ElectronCollection");
  8	HybridScCollection_ = ps.getParameter<edm::InputTag>("HybridScCollection");
  9	M5x5ScCollection_ = ps.getParameter<edm::InputTag>("M5x5ScCollection");
 10	GenMEtCollection_  = ps.getParameter<edm::InputTag>("GenMEtCollection");
 11	CaloMEtCollection_ = ps.getParameter<edm::InputTag>("CaloMEtCollection");
 12	//T1MEtCollection_ = ps.getParameter<edm::InputTag>("T1MEtCollection");
 13	PfMEtCollection_ = ps.getParameter<edm::InputTag>("PfMEtCollection");
 14	TcMEtCollection_ = ps.getParameter<edm::InputTag>("TcMEtCollection");
 15	TriggerEvent_ = ps.getParameter<edm::InputTag>("TriggerEvent");
 16	TriggerPath_ = ps.getParameter<edm::InputTag>("TriggerPath");
 17	TriggerResults_ = ps.getParameter<edm::InputTag>("TriggerResults");
 18	TriggerName_ = ps.getParameter<std::string>("TriggerName");
 19	HLTPathCheck_ = ps.getParameter<bool>("HLTPathCheck");
 20	
 21	Zevent_ = ps.getParameter<bool>("Zevent");
 22	mW_ = ps.getParameter<double>("mW");
 23        mZ_ = ps.getParameter<double>("mZ");
 24	mTPmin_ = ps.getParameter<double>("mTPmin");
 25	mTPmax_ = ps.getParameter<double>("mTPmax");
 26	BarrelEtaMax_ = ps.getParameter<double>("BarrelEtaMax");
 27	EndCapEtaMin_ = ps.getParameter<double>("EndCapEtaMin");
 28	EndCapEtaMax_ = ps.getParameter<double>("EndCapEtaMax");
 29	
 30	hyb_fCorrPSet_ = ps.getParameter<edm::ParameterSet>("hyb_fCorrPSet");
 31	m5x5_fCorrPSet_ = ps.getParameter<edm::ParameterSet>("m5x5_fCorrPSet");
 32
 33 	double CElecPtMin = ps.getParameter<double>("CElecPtMin");
 34	double CEB_siEiE = ps.getParameter<double>("CEB_sigmaIEtaIEta"); 
 35	double CEB_dPhiIn = ps.getParameter<double>("CEB_deltaPhiIn"); 
 36	double CEB_dEtaIn = ps.getParameter<double>("CEB_deltaEtaIn"); 
 37	double CEB_EcalIso = ps.getParameter<double>("CEB_EcalIso");
 38	double CEB_HcalIso = ps.getParameter<double>("CEB_HcalIso");
 39	double CEB_TrckIso = ps.getParameter<double>("CEB_TrckIso");
 40	double CEE_siEiE = ps.getParameter<double>("CEE_sigmaIEtaIEta"); 
 41	double CEE_dPhiIn = ps.getParameter<double>("CEE_deltaPhiIn"); 
 42	double CEE_dEtaIn = ps.getParameter<double>("CEE_deltaEtaIn"); 
 43	double CEE_EcalIso = ps.getParameter<double>("CEE_EcalIso");
 44	double CEE_HcalIso = ps.getParameter<double>("CEE_HcalIso");
 45	double CEE_TrckIso = ps.getParameter<double>("CEE_TrckIso");
 46
 47	CutVector_.resize(13);
 48	CutVector_[EtCut_] = CElecPtMin;
 49	CutVector_[EB_sIhIh_] = CEB_siEiE;
 50	CutVector_[EB_dPhiIn_] = CEB_dPhiIn;
 51	CutVector_[EB_dEtaIn_] = CEB_dEtaIn;
 52	CutVector_[EB_TrckIso_] = CEB_TrckIso;
 53	CutVector_[EB_EcalIso_] = CEB_EcalIso;
 54	CutVector_[EB_HcalIso_] = CEB_HcalIso;
 55	CutVector_[EE_sIhIh_] = CEE_siEiE;
 56	CutVector_[EE_dPhiIn_] = CEE_dPhiIn;
 57	CutVector_[EE_dEtaIn_] = CEE_dEtaIn;
 58	CutVector_[EE_TrckIso_] = CEE_TrckIso;
 59	CutVector_[EE_EcalIso_] = CEE_EcalIso;
 60	CutVector_[EE_HcalIso_] = CEE_HcalIso;
 61	
 62	for(std::vector<double>::const_iterator it = CutVector_.begin(); it != CutVector_.end(); ++it)
 63	{
 64		edm::LogDebug_("","",101)<<"CutVector_ = "<< *it;
 65	}
 66}
 67
 68ErsatzMEt::~ErsatzMEt()
 69{
 70}
 71
 72// ------------ method called once each job just before starting event loop  ------------
 73void ErsatzMEt::beginJob()
 74{
 75	edm::Service<TFileService> fs;
 76
 77	t_ = fs->make<TTree>("ErsatzMEt", "Data on ErsatzMEt");
 78	
 79	edm::LogDebug_("","", 75)<<"Creating Ersatz MEt branches.";
 80	
 81	t_->Branch("nTags", &nTags_, "nTags/I");
 82	t_->Branch("nProbes", &nProbes_, "nProbes/I");
 83	
 84	t_->Branch("ErsatzV1CaloMEt", ErsatzV1CaloMEt_, "ErsatzV1CaloMEt[4]/D");
 85	t_->Branch("ErsatzV1CaloMt", ErsatzV1CaloMt_, "ErsatzV1CaloMt[4]/D");
 86	t_->Branch("ErsatzV1CaloMEtPhi", ErsatzV1CaloMEtPhi_, "ErsatzV1CaloMEtPhi[4]/D");
 87	t_->Branch("ErsatzV2CaloMEt", ErsatzV2CaloMEt_, "ErsatzV2CaloMEt[4]/D");
 88	t_->Branch("ErsatzV2CaloMEtPhi", ErsatzV2CaloMEtPhi_, "ErsatzV2CaloMEtPhi[4]/D");
 89	t_->Branch("ErsatzV2CaloMt", ErsatzV2CaloMt_, "ErsatzV2CaloMt[4]/D");
 90	t_->Branch("ErsatzV3CaloMEt", ErsatzV3CaloMEt_, "ErsatzV3CaloMEt[4]/D");
 91	t_->Branch("ErsatzV3CaloMEtPhi", ErsatzV3CaloMEtPhi_, "ErsatzV3CaloMEtPhi[4]/D");
 92	t_->Branch("ErsatzV3CaloMt", ErsatzV3CaloMt_, "ErsatzV3CaloMt[4]/D");
 93	t_->Branch("ErsatzV4CaloMEt", ErsatzV4CaloMEt_, "ErsatzV4CaloMEt[4]/D");
 94	t_->Branch("ErsatzV4CaloMEtPhi", ErsatzV4CaloMEtPhi_, "ErsatzV4CaloMEtPhi[4]/D");
 95	t_->Branch("ErsatzV4CaloMt", ErsatzV4CaloMt_, "ErsatzV4CaloMt[4]/D");
 96	t_->Branch("ErsatzV1T1MEt", ErsatzV1T1MEt_, "ErsatzV1T1MEt[4]/D");
 97	t_->Branch("ErsatzV1T1Mt", ErsatzV1T1Mt_, "ErsatzV1T1Mt[4]/D");
 98	t_->Branch("ErsatzV1T1MEtPhi", ErsatzV1T1MEtPhi_, "ErsatzV1T1MEtPhi[4]/D");
 99	t_->Branch("ErsatzV1PfMEt", ErsatzV1PfMEt_, "ErsatzV1PfMEt[4]/D");
100	t_->Branch("ErsatzV1PfMt", ErsatzV1PfMt_, "ErsatzV1PfMt[4]/D");
101	t_->Branch("ErsatzV1PfMEtPhi", ErsatzV1TcMEtPhi_, "ErsatzV1PfMEtPhi[4]/D");
102	t_->Branch("ErsatzV1TcMEt", ErsatzV1TcMEt_, "ErsatzV1TcMEt[4]/D");
103	t_->Branch("ErsatzV1TcMt", ErsatzV1TcMt_, "ErsatzV1TcMt[4]/D");
104	t_->Branch("ErsatzV1TcMEtPhi", ErsatzV1TcMEtPhi_, "ErsatzV1TcMEtPhi[4]/D");
105
106	t_->Branch("CaloMEt", &CaloMEt_, "CaloMEt/D");
107	t_->Branch("CaloMEtphi", &CaloMEtphi_, "CaloMEtphi/D");
108	t_->Branch("T1MEt", &T1MEt_, "T1MEt/D");
109	t_->Branch("T1MEtphi", &T1MEtphi_, "T1MEtphi/D");
110	t_->Branch("PfMEt", &PfMEt_, "PfMEt/D");
111	t_->Branch("PfMEtphi", &PfMEtphi_, "PfMEtphi/D");
112	t_->Branch("TcMEt", &TcMEt_, "TcMEt/D");
113	t_->Branch("TcMEtphi", &TcMEtphi_, "TcMEtphi/D");
114
115	edm::LogDebug_("","", 91)<<"Creating electron branches.";
116	t_->Branch("tag_q", tag_q_,"tag_q[4]/I");
117	t_->Branch("tag_pt", tag_pt_,"tag_pt[4]/D");
118	t_->Branch("tag_eta", tag_eta_,"tag_eta[4]/D");
119	t_->Branch("tag_phi", tag_phi_,"tag_phi[4]/D");
120	t_->Branch("tag_sIhIh", tag_sIhIh_, "tag_sIhIh[4]/D");
121	t_->Branch("tag_dPhiIn", tag_dPhiIn_, "tag_dPhiIn[4]/D");
122	t_->Branch("tag_dEtaIn", tag_dEtaIn_, "tag_dEtaIn[4]/D");
123	t_->Branch("tag_trckIso", tag_trckIso_,"tag_trckIso[4]/D");
124	t_->Branch("tag_ecalIso", tag_ecalIso_,"tag_ecalIso[4]/D");
125	t_->Branch("tag_hcalIso", tag_hcalIso_,"tag_hcalIso[4]/D");
126	t_->Branch("tag_e2x5Max", tag_e2x5Max_,"tag_e2x5Max[4]/D");
127	t_->Branch("tag_e1x5Max", tag_e1x5Max_,"tag_e1x5Max[4]/D");
128	t_->Branch("tag_e5x5", tag_e5x5_,"tag_e5x5[4]/D");
129	t_->Branch("tag_hoe", tag_hoe_,"tag_hoe[4]/D");
130	t_->Branch("tag_eop", tag_eop_,"tag_eop[4]/D");
131	t_->Branch("tag_pin", tag_pin_,"tag_pin[4]/D");
132	t_->Branch("tag_pout", tag_pout_,"tag_pout[4]/D");
133	t_->Branch("tag_rescPt", tag_rescPt_, "tag_rescPt[4]/D");
134	t_->Branch("tag_rescEta", tag_rescEta_, "tag_rescEta[4]/D");
135	t_->Branch("tag_rescPhi", tag_rescPhi_, "tag_rescPhi[4]/D");
136	
137	edm::LogDebug_("","", 103)<<"Creating ersatz neutrino branches.";
138	t_->Branch("probe_q", probe_q_,"probe_q[4]/I");
139	t_->Branch("probe_pt", probe_pt_,"probe_pt[4]/D");
140	t_->Branch("probe_eta", probe_eta_,"probe_eta[4]/D");
141	t_->Branch("probe_phi", probe_phi_,"probe_phi[4]/D");
142	t_->Branch("probe_sIhIh", probe_sIhIh_, "probe_sIhIh[4]/D");
143	t_->Branch("probe_dPhiIn", probe_dPhiIn_, "probe_dPhiIn[4]/D");
144	t_->Branch("probe_dEtaIn", probe_dEtaIn_, "probe_dEtaIn[4]/D");
145	t_->Branch("probe_trckIso", probe_trckIso_,"probe_trckIso[4]/D");
146	t_->Branch("probe_ecalIso", probe_ecalIso_,"probe_ecalIso[4]/D");
147	t_->Branch("probe_hcalIso", probe_hcalIso_,"probe_hcalIso[4]/D");
148	t_->Branch("probe_e2x5Max", probe_e2x5Max_,"probe_e2x5Max[4]/D");
149	t_->Branch("probe_e1x5Max", probe_e1x5Max_,"probe_e1x5Max[4]/D");
150	t_->Branch("probe_e5x5", probe_e5x5_,"probe_e5x5[4]/D");
151	t_->Branch("probe_hoe", probe_hoe_,"probe_hoe[4]/D");
152	t_->Branch("probe_eop", probe_eop_,"probe_eop[4]/D");
153	t_->Branch("probe_pin", probe_pin_,"probe_pin[4]/D");
154	t_->Branch("probe_pout", probe_pout_,"probe_pout[4]/D");
155	t_->Branch("probe_rescPt", probe_rescPt_, "probe_rescPt[4]/D");
156	t_->Branch("probe_rescEta", probe_rescEta_, "probe_rescEta[4]/D");
157	t_->Branch("probe_rescPhi", probe_rescPhi_, "probe_rescPhi[4]/D");
158
159	t_->Branch("Z_m", Z_m_, "Z_m[4]/D");
160	t_->Branch("Z_pt", Z_pt_, "Z_pt[4]/D");
161	t_->Branch("Z_eta", Z_eta_, "Z_eta[4]/D");
162	t_->Branch("Z_y", Z_y_, "Z_y[4]/D");
163	t_->Branch("Z_phi", Z_phi_, "Z_phi[4]/D");
164	t_->Branch("Z_rescM", Z_rescM_, "Z_rescM[4]/D");
165	t_->Branch("Z_rescPt", Z_rescPt_, "Z_rescPt[4]/D");
166	t_->Branch("Z_rescEta", Z_rescEta_, "Z_rescEta[4]/D");
167	t_->Branch("Z_rescY", Z_rescY_, "Z_rescY[4]/D");
168	t_->Branch("Z_rescPhi", Z_rescPhi_, "Z_rescPhi[4]/D");
169	t_->Branch("Z_probe_dPhi",Z_probe_dPhi_,"Z_probe_dPhi[4]/D"); 
170
171	t_->Branch("probe_sc_pt",probe_sc_pt_, "probe_sc_pt[4]/D");
172	t_->Branch("probe_sc_eta",probe_sc_eta_, "probe_sc_eta[4]/D");
173	t_->Branch("probe_sc_phi", probe_sc_phi_, "probe_sc_phi[4]/D");
174	t_->Branch("probe_sc_E",probe_sc_E_, "probe_sc_E[4]/D");
175	t_->Branch("probe_sc_rawE",probe_sc_rawE_, "probe_sc_rawE[4]/D");
176	t_->Branch("probe_sc_nClus", probe_sc_nClus_, "probe_sc_nClus[4]/D");
177	t_->Branch("probe_scV2_E",probe_scV2_E_, "probe_scV2_E[4]/D");
178	t_->Branch("probe_scV3_E",probe_scV3_E_, "probe_scV3_E[4]/D");
179	t_->Branch("probe_scV4_E",probe_scV4_E_, "probe_scV4_E[4]/D");
180	t_->Branch("probe_d_MCE_SCE", probe_d_MCE_SCE_, "probe_d_MCE_SCE[4]/D");
181
182	t_->Branch("ErsatzV1_Mesc", ErsatzV1_Mesc_, "ErsatzV1_Mesc[4]/D");
183	t_->Branch("ErsatzV1_rescMesc", ErsatzV1_rescMesc_, "ErsatzV1_rescMesc[4]/D");
184
185	t_->Branch("McElec_nFinal", &McElec_nFinal_, "McElec_nFinal/I");
186
187	if(Zevent_){
188		t_->Branch("McZ_m", &McZ_m_, "McZ_m/D");
189		t_->Branch("McZ_rescM", &McZ_rescM_, "McZ_rescM/D");
190		t_->Branch("McZ_Pt", &McZ_pt_, "McZ_Pt/D");
191		t_->Branch("McZ_y", &McZ_y_, "McZ_y/D");
192		t_->Branch("McZ_Eta", &McZ_eta_, "McZ_Eta/D");
193		t_->Branch("McZ_Phi", &McZ_phi_, "McZ_Phi/D");
194		t_->Branch("McZ_rescPt", &McZ_rescPt_, "McZ_Pt/D");
195		t_->Branch("McZ_rescY", &McZ_rescY_, "McZ_rescY/D");
196		t_->Branch("McZ_rescEta", &McZ_rescEta_, "McZ_Eta/D");
197		t_->Branch("McZ_rescPhi", &McZ_rescPhi_, "McZ_Phi/D");
198		t_->Branch("McElec_nZmum", &McElec_nZmum_, "McElec_nZmum/I");
199		t_->Branch("McElec_eta", &McElec_eta_, "McElec_eta[4]/D");
200		t_->Branch("McElec_pt", &McElec_pt_, "McElec_pt[4]/D");
201		t_->Branch("McElec_phi", &McElec_phi_, "McElec_phi[4]/D");
202		t_->Branch("McElec_rescEta", &McElec_rescEta_, "McElec_rescEta[4]/D");
203		t_->Branch("McElec_rescPhi", &McElec_rescPhi_, "McElec_rescPhi[4]/D");
204		t_->Branch("McElec_rescPt", &McElec_rescPt_, "McElec_rescPt[4]/D");
205		t_->Branch("McProbe_eta", &McProbe_eta_, "McProbe_eta[4]/D");
206		t_->Branch("McProbe_pt", &McProbe_pt_, "McProbe_pt[4]/D");
207		t_->Branch("McProbe_phi", &McProbe_phi_, "McProbe_phi[4]/D");
208		t_->Branch("McProbe_rescEta", &McProbe_rescEta_, "McProbe_rescEta[4]/D");
209		t_->Branch("McProbe_rescPt", &McProbe_rescPt_, "McProbe_rescPt[4]/D");
210		t_->Branch("McProbe_rescPhi", &McProbe_rescPhi_, "McProbe_rescPhi[4]/D");
211		t_->Branch("McElecProbe_dPhi", &McElecProbe_dPhi_, "McElecProbe_dPhi/D");
212		t_->Branch("McElecProbe_dEta", &McElecProbe_dEta_, "McElecProbe_dEta/D");
213		t_->Branch("McElecProbe_dR", &McElecProbe_dR_, "McElecProbe_dR/D");
214	}
215
216}
217
218void ErsatzMEt::analyze(const edm::Event& evt, const edm::EventSetup& es)
219{
220
221       es.get<CaloGeometryRecord>().get(geoHandle_);
222       es.get<CaloTopologyRecord>().get(pTopology_);
223
224	edm::LogDebug_("","", 151)<<"Initialising variables.";
225	nTags_ = -99; nProbes_ = -99;
226	CaloMEt_ = -99.; CaloMEtphi_ = -99.;
227	T1MEt_ = -99.; T1MEtphi_ = -99.;
228	PfMEt_ = -99.; PfMEtphi_ = -99.;
229	TcMEt_ = -99.; TcMEtphi_ = -99.;	
230	if(Zevent_)
231	{
232		McZ_m_ = -99.; McZ_pt_ = -99.; McZ_y_ = -99.; McZ_eta_ = -99.; McZ_phi_ = -99.;
233		McZ_rescM_ = -99.; McZ_rescPt_ = -99.; McZ_rescY_ = -99.; McZ_rescEta_ = -99.; McZ_rescPhi_ = -99.;
234		McElec_nZmum_ = -99; McElec_nFinal_ = -99;
235	}
236
237	for(int i = 0; i < nEntries_arr_; ++i)
238	{
239		tag_q_[i] = -99;
240		tag_pt_[i] = -99.; tag_eta_[i] = -99.; tag_phi_[i] = -99.;
241		tag_rescPt_[i] = -99.; tag_rescEta_[i] = -99.; tag_rescPhi_[i] = -99.;
242		tag_trckIso_[i] = -99.; tag_ecalIso_[i] = -99.; tag_hcalIso_[i] = -99.;
243		tag_sIhIh_[i] = -99.; tag_dPhiIn_[i] = -99.; tag_dEtaIn_[i] = -99.;
244		tag_e5x5_[i] = -99.; tag_e2x5Max_[i] = -99.; tag_e1x5Max_[i] = -99.;
245		tag_hoe_[i] = -99.; tag_eop_[i] = -99.; tag_pin_[i] = -99.; tag_pout_[i] = -99.;
246
247		probe_q_[i] = -99;
248		probe_pt_[i] = -99.; probe_eta_[i] = -99.; probe_phi_[i] = -99.;
249		probe_rescPt_[i] = -99.; probe_rescEta_[i] = -99.; probe_rescPhi_[i] = -99.;
250		probe_trckIso_[i] = -99.; probe_ecalIso_[i] = -99.; probe_hcalIso_[i] = -99.;
251		probe_sIhIh_[i] = -99.; probe_dPhiIn_[i] = -99.; probe_dEtaIn_[i] = -99.;
252		probe_e5x5_[i] = -99.; probe_e2x5Max_[i] = -99.; probe_e1x5Max_[i] = -99.;
253		probe_hoe_[i] = -99.; probe_eop_[i] = -99.; probe_pin_[i] = -99.; probe_pout_[i] = -99.;
254		
255		Z_pt_[i] = -99.; Z_y_[i] = -99.; Z_eta_[i] = -99.; Z_phi_[i] = -99.; Z_m_[i] = -99.;
256		Z_rescPt_[i] = -99.; Z_rescY_[i] = -99.; Z_rescEta_[i] = -99.; Z_rescPhi_[i] = -99.; Z_rescM_[i] = -99.; Z_probe_dPhi_[i] = -99.;
257		
258		ErsatzV1_Mesc_[i] = -99.; ErsatzV1_rescMesc_[i] = -99.;
259		ErsatzV2_Mesc_[i] = -99.; ErsatzV2_rescMesc_[i] = -99.;
260		ErsatzV3_Mesc_[i] = -99.; ErsatzV3_rescMesc_[i] = -99.;
261		ErsatzV4_Mesc_[i] = -99.; ErsatzV4_rescMesc_[i] = -99.;
262		ErsatzV1CaloMEt_[i] = -99.; ErsatzV1CaloMt_[i] = -99.; ErsatzV1CaloMEtPhi_[i] = -99.; 
263		ErsatzV2CaloMEt_[i] = -99.; ErsatzV2CaloMt_[i] = -99.; ErsatzV2CaloMEtPhi_[i] = -99.; 
264		ErsatzV3CaloMEt_[i] = -99.; ErsatzV3CaloMt_[i] = -99.; ErsatzV3CaloMEtPhi_[i] = -99.; 
265		ErsatzV4CaloMEt_[i] = -99.; ErsatzV4CaloMt_[i] = -99.; ErsatzV4CaloMEtPhi_[i] = -99.; 
266		ErsatzV1T1MEt_[i] = -99.; ErsatzV1T1Mt_[i] = -99.; ErsatzV1T1MEtPhi_[i] = -99.; 
267		ErsatzV1PfMEt_[i] = -99.; ErsatzV1PfMt_[i] = -99.; ErsatzV1PfMEtPhi_[i] = -99.; 
268		ErsatzV1TcMEt_[i] = -99.; ErsatzV1TcMt_[i] = -99.; ErsatzV1TcMEtPhi_[i] = -99.; 
269		
270		probe_sc_pt_[i] = -99.; probe_sc_eta_[i] = -99.; probe_sc_phi_[i] = -99.;
271		probe_sc_E_[i] = -99.; probe_sc_rawE_[i] = -99.;
272		probe_scV2_E_[i] = -99.;
273		probe_scV3_E_[i] = -99.;
274		probe_scV4_E_[i] = -99.;
275
276		if(Zevent_)
277		{
278			McElec_pt_[i] = -99.; McElec_eta_[i] = -99.; McElec_phi_[i] = -99.;
279			McElec_rescPt_[i] = -99.; McElec_rescEta_[i] = -99.; McElec_rescPhi_[i] = -99.;
280			McProbe_pt_[i] = -99.; McProbe_eta_[i] = -99.; McProbe_phi_[i] = -99.;
281			McProbe_rescPt_[i] = -99.; McProbe_rescEta_[i] = -99.; McProbe_rescPhi_[i] = -99.;
282			McElecProbe_dPhi_[i] = -99.; McElecProbe_dEta_[i] = -99.; McElecProbe_dR_[i] = -99.;
283		}
284
285		edm::LogDebug_("","",180)<<"Initialisation of array index "<< i <<" completed.";
286	}	
287	//Get Collections
288	edm::Handle<reco::GenParticleCollection> pGenPart;
289	try
290	{
291		evt.getByLabel(MCTruthCollection_, pGenPart);
292	}catch(cms::Exception &ex)
293	{	
294		edm::LogError("analyze") <<"Can't get collection with label "<< MCTruthCollection_.label();
295	} 
296	edm::Handle<reco::GsfElectronCollection> pElectrons;
297	try
298	{
299		evt.getByLabel(ElectronCollection_, pElectrons);
300	}catch(cms::Exception &ex)
301	{
302		edm::LogError("analyze") <<"Can't get collection with label "<< ElectronCollection_.label();
303	}
304	edm::Handle<reco::SuperClusterCollection> pHybrid;
305	try
306	{
307		evt.getByLabel(HybridScCollection_, pHybrid);
308	}catch(cms::Exception &ex)
309	{
310		edm::LogError("analyze") <<"Can't get collection with label "<< HybridScCollection_.label();
311	}
312	edm::Handle<reco::SuperClusterCollection> pM5x5;
313	try
314	{
315		evt.getByLabel(M5x5ScCollection_, pM5x5);
316	}catch(cms::Exception &ex)
317	{
318		edm::LogError("analyze") <<"Can't get collection with label "<< M5x5ScCollection_.label();
319	}
320	edm::Handle<reco::CaloMETCollection> pCaloMEt;
321	try
322	{
323		evt.getByLabel(CaloMEtCollection_, pCaloMEt);
324	}catch(cms::Exception &ex)
325	{
326		edm::LogError("analyze")<<"Can't get collection with label "<< CaloMEtCollection_.label();
327	}
328	edm::Handle<reco::METCollection> pT1MEt;
329//	try
330//	{
331//		evt.getByLabel(T1MEtCollection_, pT1MEt);
332//	}catch(cms::Exception &ex)
333//	{
334//		edm::LogError("analyze")<<"Can't get collection with label "<< T1MEtCollection_.label();
335//	}
336	edm::Handle<reco::PFMETCollection> pPfMEt;
337	try
338	{
339		evt.getByLabel(PfMEtCollection_, pPfMEt);
340	}catch(cms::Exception &ex)
341	{
342		edm::LogError("analyze")<<"Can't get collection with label "<< PfMEtCollection_.label();
343	}
344	edm::Handle<reco::METCollection> pTcMEt;
345	try
346	{
347		evt.getByLabel(TcMEtCollection_, pTcMEt);
348	}catch(cms::Exception &ex)
349	{
350		edm::LogError("analyze")<<"Can't get collection with label "<< TcMEtCollection_.label();
351	}
352	edm::Handle<reco::GenMETCollection> pGenMEt;
353	try
354	{
355		evt.getByLabel(GenMEtCollection_, pGenMEt);
356	}catch(cms::Exception &ex)
357	{
358		edm::LogError("analyze")<<"Can't get collection with label "<< GenMEtCollection_.label();
359	}
360	edm::Handle<edm::TriggerResults> pTriggerResults;
361	try
362	{
363		evt.getByLabel(TriggerResults_, pTriggerResults);
364	}catch(cms::Exception &ex)
365	{
366		edm::LogError("analyze")<<"Cant get collection with label "<< TriggerResults_.label();
367	}
368	edm::Handle<trigger::TriggerEvent> pHLT;
369	try
370	{
371		evt.getByLabel(TriggerEvent_, pHLT);
372	}catch(cms::Exception &ex)
373	{
374		edm::LogError("analyze")<<"Can't get collection with label "<< TriggerEvent_.label();
375	}
376
377	std::vector<math::XYZTLorentzVector>McElecs,McElecsFinalState;
378	std::vector<math::XYZTLorentzVector> McElecsResc;
379	if(Zevent_)
380	{
381		edm::LogDebug_("","",289)<<"Analysing MC properties.";
382		const reco::GenParticleCollection *McCand = pGenPart.product();
383		math::XYZTLorentzVector Zboson, RescZboson, McElec1, McElec2; 
384		for(reco::GenParticleCollection::const_iterator McP = McCand->begin(); McP != McCand->end(); ++McP)
385		{
386			const reco::Candidate* mum = McP->mother();
387			if(abs(McP->pdgId())==11 && abs(mum->pdgId()) == 23)
388			{
389				McElecs.push_back(McP->p4());
390				if(abs(mum->pdgId() == 23)) Zboson = mum->p4();
391				
392				std::cout <<"Found electron, ID = "<< McP->pdgId() <<"\t status = "<< McP->status()<<std::endl;
393				if(McP->status() != 1)
394				{
395					const reco::Candidate* McPD = McP->daughter(0);
396					McPD = McPD->mother();
397					while(McPD->status() != 1)
398					{
399						int n = McPD->numberOfDaughters();
400						std::cout<< McPD->pdgId() << " : status = "<<McPD->status()
401							<<"\tNumber of Daughters = "<< n <<std::endl;
402						for(int j = 0; j < n; ++ j) 
403						{
404							const reco::Candidate *d = McPD->daughter( j );
405							std::cout <<"Daughter "<< j <<"\t id = "<< d->pdgId() << std::endl;
406							if(abs(d->pdgId()) == 11)
407							{ 
408								McPD = d;
409								break;
410							}
411						}
412					}
413					std::cout<< McPD->pdgId() << " : status = "<<McPD->status()<<"\tAdding to vector!"<<std::endl;
414					McElecsFinalState.push_back(McPD->p4());
415				}else McElecsFinalState.push_back(McP->p4());
416			}
417		}
418		McZ_m_ = Zboson.M(); McZ_pt_ = Zboson.Pt(); McZ_phi_ = Zboson.Phi(); McZ_eta_ = Zboson.Eta(); McZ_y_ = Zboson.Rapidity();
419		McElec_nZmum_ =McElecs.size();
420		McElec_nFinal_ =McElecsFinalState.size();
421		edm::LogDebug_("","",309)<<"MC electrons with Z mother = "<< McElec_nZmum_ 
422						<<"\tFinal state MC electrons = "<< McElec_nFinal_;
423		
424		McElecsResc.resize(2);
425//		RescZboson.SetCoordinates(Zboson.Px(), Zboson.Py(), Zboson.Pz(), sqrt(Zboson.P2()+(mW_*mW_*Zboson.M2())/(mZ_*mZ_)));
426		RescZboson.SetCoordinates(Zboson.Px()*mW_/mZ_, Zboson.Py()*mW_/mZ_, Zboson.Pz()*mW_/mZ_, Zboson.E()*mW_/mZ_);
427		McZ_rescM_ = RescZboson.M(); McZ_rescPt_ = RescZboson.Pt(); McZ_rescEta_ = RescZboson.Eta(); McZ_rescPhi_ = RescZboson.Phi();
428		McZ_rescY_ = RescZboson.Rapidity();
429		ROOT::Math::Boost CoMBoost(Zboson.BoostToCM());
430		
431		math::XYZTLorentzVector RescMcElec0 = CoMBoost(McElecsFinalState[0]);
432		math::XYZTLorentzVector RescMcElec1 = CoMBoost(McElecsFinalState[1]);
433		RescMcElec0 *= mW_/mZ_;
434		RescMcElec1 *= mW_/mZ_;
435
436		double E_W = RescZboson.E();
437		ROOT::Math::Boost BackToLab(RescZboson.Px()/E_W, RescZboson.Py()/E_W, RescZboson.Pz()/E_W);
438
439		RescMcElec0 = BackToLab(RescMcElec0);
440//		RndmMcElec_Rescaled_pt_ = RescMcElec0.Pt();
441//		RndmMcElec_Rescaled_eta_ = RescMcElec0.Eta();
442//		RndmMcElec_Rescaled_phi_ = RescMcElec0.Phi();
443	
444		RescMcElec1 = BackToLab(RescMcElec1);
445//		OthrMcElec_Rescaled_pt_ = RescMcElec1.Pt();
446//		OthrMcElec_Rescaled_eta_ = RescMcElec1.Eta();
447//		OthrMcElec_Rescaled_phi_ = RescMcElec1.Phi();
448		McElecsResc[0] = RescMcElec0;
449		McElecsResc[1] = RescMcElec1;
450		math::XYZTLorentzVector sum = RescMcElec1+RescMcElec0;
451		edm::LogDebug_("","", 307)<<"McElecsResc[0] + McElecsResc[1] = ("<<sum.Px()<<", "<<sum.Py()<<", "
452						<<sum.Pz()<<", "<<sum.E()<<")";
453	}	
454
455	const edm::TriggerResults* HltRes = pTriggerResults.product();
456	const edm::TriggerNames & triggerNames = evt.triggerNames(*HltRes);
457	if(HLTPathCheck_)
458	{
459		for(unsigned int itrig = 0; itrig < HltRes->size(); ++itrig)
460		{
461			std::string nom = triggerNames.triggerName(itrig);
462			edm::LogInfo("")<< itrig <<" : Name = "<< nom <<"\t Accepted = "<< HltRes->accept(itrig);
463		}
464	}
465	if(HltRes->accept(34) ==0) edm::LogError("")<<"Event did not pass "<< triggerNames.triggerName(34)<<"!";
466	if(HltRes->accept(34) !=0)
467	{
468	std::vector<reco::GsfElectronRef> UniqueElectrons;
469	UniqueElectrons = uniqueElectronFinder(pElectrons);
470	edm::LogDebug_("","ErsatzMEt",192)<<"Unique electron size = "<<UniqueElectrons.size();
471	std::vector<reco::GsfElectronRef> SelectedElectrons;
472	const unsigned int fId = pHLT->filterIndex(TriggerPath_);
473	std::cout << "Filter Id = " << fId << std::endl;
474	SelectedElectrons = electronSelector(UniqueElectrons, pHLT, fId, CutVector_);
475	nTags_ = SelectedElectrons.size();
476	edm::LogDebug_("","ErsatzMEt",197)<<"Selected electron size = "<<nTags_;
477
478	iComb_ = 0;
479	if(Zevent_)
480	{
481	//Match MC electrons to the selected electrons and store some of their properties in the tree.
482	//The properties of the other MC electron (i.e. that not selected) are also stored.
483		for(std::vector<reco::GsfElectronRef>::const_iterator elec = SelectedElectrons.begin(); 
484								elec != SelectedElectrons.end(); ++elec)
485		{
486			for(int m = 0; m < 2; ++m)
487			{
488				double dRLimit = 99.;
489				double dR = reco::deltaR(McElecs[m], *(*elec)); 
490				if(dR < dRLimit)
491				{
492					dRLimit = dR;
493					McElec_pt_[iComb_] = McElecs[m].pt();
494					McElec_eta_[iComb_] = McElecs[m].eta();
495					McElec_rescPt_[iComb_] = McElecsResc[m].pt();
496					McElec_rescEta_[iComb_] = McElecsResc[m].eta();
497				} 
498			}
499		}
500	}
501
502	std::map<reco::GsfElectronRef, reco::GsfElectronRef> TagProbePairs;
503	TagProbePairs = probeFinder(SelectedElectrons, pElectrons);
504	nProbes_ = TagProbePairs.size();	
505	edm::LogDebug_("", "ErsatzMEt", 209)<<"Number of tag-probe pairs = "<< TagProbePairs.size();
506
507	if(!TagProbePairs.empty())
508	{
509		const reco::CaloMETCollection* caloMEtCollection = pCaloMEt.product();
510        	const reco::MET calomet = *(caloMEtCollection->begin());
511		CaloMEt_ = calomet.pt();
512		CaloMEtphi_ = calomet.phi();
513
514		//const reco::METCollection* t1MEtCollection = pT1MEt.product();
515        	//const reco::MET t1met = *(t1MEtCollection->begin());
516		//T1MEt_ = t1met.pt();
517		//T1MEtphi_ = t1met.phi();
518		
519		const reco::PFMETCollection* pfMEtCollection = pPfMEt.product();
520        	const reco::MET pfmet = *(pfMEtCollection->begin());
521		PfMEt_ = pfmet.pt();
522		PfMEtphi_ = pfmet.phi();
523		
524		const reco::METCollection* tcMEtCollection = pTcMEt.product();
525        	const reco::MET tcmet = *(tcMEtCollection->begin());
526		TcMEt_ = tcmet.pt();
527		TcMEtphi_ = tcmet.phi();
528		
529		reco::MET ersatzMEt;
530
531		for(std::map<reco::GsfElectronRef, reco::GsfElectronRef>::const_iterator it = TagProbePairs.begin();
532			it != TagProbePairs.end(); ++it)
533		{
534			edm::LogDebug_("","DelendumLoop", 293)<<"iComb_ = "<< iComb_;
535			tag_q_[iComb_] = it->first->charge();
536			edm::LogDebug_("","",360)<<"tag charge = "<< tag_q_[iComb_];
537			tag_pt_[iComb_] = it->first->pt();
538			tag_eta_[iComb_] = it->first->eta();
539			tag_phi_[iComb_] = it->first->phi();
540			edm::LogDebug_("","ErsatzMEt", 364)<<"tag pt = "<< tag_pt_[iComb_] 
541					<<"\teta = "<< tag_eta_[iComb_]<<"\tphi = "<< tag_phi_[iComb_];
542                        tag_trckIso_[iComb_] = it->first->isolationVariables03().tkSumPt;
543                        tag_ecalIso_[iComb_] = it->first->isolationVariables04().ecalRecHitSumEt;
544                        tag_hcalIso_[iComb_] = it->first->isolationVariables04().hcalDepth1TowerSumEt 
545						+ it->first->isolationVariables04().hcalDepth2TowerSumEt;
546			edm::LogDebug_("","ErsatzMEt", 370)<<"tag trackiso = "<< tag_trckIso_[iComb_] 
547					<<"\tecaliso = "<< tag_ecalIso_[iComb_]<<"\thcaliso = "<< tag_hcalIso_[iComb_];
548			tag_sIhIh_[iComb_] = it->first->scSigmaIEtaIEta();
549			tag_dPhiIn_[iComb_] = it->first->deltaPhiSuperClusterTrackAtVtx();
550			tag_dEtaIn_[iComb_] = it->first->deltaEtaSuperClusterTrackAtVtx();
551			edm::LogDebug_("","ErsatzMEt", 245)<<"tag sIhIh = "<< tag_sIhIh_[iComb_] 
552					<<"\tdPhiIn = "<< tag_dPhiIn_[iComb_]<<"\tdEtaIn = "<< tag_dEtaIn_[iComb_];
553			tag_e5x5_[iComb_] = it->first->scE5x5();
554			tag_e2x5Max_[iComb_] = it->first->scE2x5Max();
555			tag_e2x5Max_[iComb_] = it->first->scE1x5();
556			edm::LogDebug_("","ErsatzMEt", 245)<<"tag e5x5 = "<< tag_e5x5_[iComb_] 
557					<<"\te2x5Max = "<< tag_e2x5Max_[iComb_]<<"\te1x5Max = "<< tag_e1x5Max_[iComb_];
558			tag_hoe_[iComb_] = it->first->hadronicOverEm();
559			tag_eop_[iComb_] = it->first->eSuperClusterOverP();
560			tag_pin_[iComb_] = it->first->trackMomentumAtVtx().R();
561			tag_pout_[iComb_] = it->first->trackMomentumOut().R();
562			edm::LogDebug_("","ErsatzMEt", 245)<<"tag hoe = "<<tag_hoe_[iComb_]<<"\tpoe = "<<tag_eop_[iComb_] 
563					<<"\tpin = "<< tag_pin_[iComb_]<<"\tpout = "<< tag_pout_[iComb_];
564			probe_q_[iComb_] = it->first->charge();
565			edm::LogDebug_("","",360)<<"probe charge = "<< probe_q_[iComb_];
566			probe_pt_[iComb_] = it->second->pt();
567			probe_eta_[iComb_] = it->second->eta();
568			probe_phi_[iComb_] = it->second->phi();
569			edm::LogDebug_("","ErsatzMEt", 245)<<"probe pt = "<< probe_pt_[iComb_] 
570					<<"\teta = "<< probe_eta_[iComb_]<<"\tphi = "<< probe_phi_[iComb_];
571                        probe_trckIso_[iComb_] = it->second->isolationVariables03().tkSumPt;
572                        probe_ecalIso_[iComb_] = it->second->isolationVariables04().ecalRecHitSumEt;
573                        probe_hcalIso_[iComb_] = it->second->isolationVariables04().hcalDepth1TowerSumEt 
574						+ it->second->isolationVariables04().hcalDepth2TowerSumEt;
575			edm::LogDebug_("","ErsatzMEt", 245)<<"probe trackiso = "<< probe_trckIso_[iComb_] 
576					<<"\tecaliso = "<< probe_ecalIso_[iComb_]<<"\thcaliso = "<< probe_phi_[iComb_];
577			probe_sIhIh_[iComb_] = it->second->scSigmaIEtaIEta();
578			probe_dPhiIn_[iComb_] = it->second->deltaPhiSuperClusterTrackAtVtx();
579			probe_dEtaIn_[iComb_] = it->second->deltaEtaSuperClusterTrackAtVtx();
580			edm::LogDebug_("","ErsatzMEt", 245)<<"probe sIhIh = "<< probe_sIhIh_[iComb_] 
581					<<"\tdPhiIn = "<< probe_dPhiIn_[iComb_]<<"\tdEtaIn = "<< probe_dEtaIn_[iComb_];
582			probe_e5x5_[iComb_] = it->second->scE5x5();
583			probe_e2x5Max_[iComb_] = it->second->scE2x5Max();
584			probe_e2x5Max_[iComb_] = it->second->scE1x5();
585			edm::LogDebug_("","ErsatzMEt", 245)<<"probe e5x5 = "<< probe_e5x5_[iComb_] 
586					<<"\te2x5Max = "<< probe_e2x5Max_[iComb_]<<"\te1x5Max = "<< probe_e1x5Max_[iComb_];
587			probe_hoe_[iComb_] = it->second->hadronicOverEm();
588			probe_eop_[iComb_] = it->second->eSuperClusterOverP();
589			probe_pin_[iComb_] = it->second->trackMomentumAtVtx().R();
590			probe_pout_[iComb_] = it->second->trackMomentumOut().R();
591			edm::LogDebug_("","ErsatzMEt", 245)<<"probe hoe = "<<probe_hoe_[iComb_]<<"\tpoe = "<<probe_eop_[iComb_] 
592					<<"\tpin = "<< probe_pin_[iComb_]<<"\tpout = "<< probe_pout_[iComb_];
593			
594			double dRLimit = 0.2;
595			for(unsigned int mcEId = 0; mcEId < McElecs.size(); ++mcEId)
596			{
597//				double dR = reco::deltaR((*(*mcEl)), probeVec); 
598				double dR = reco::deltaR(McElecs[mcEId], it->second->p4()); 
599				if(dR < dRLimit)
600				{
601					dRLimit = dR;
602					McProbe_pt_[iComb_] = McElecs[mcEId].pt();
603					McProbe_eta_[iComb_] = McElecs[mcEId].eta();
604					McProbe_phi_[iComb_] = McElecs[mcEId].phi();
605					McProbe_rescPt_[iComb_] = McElecsResc[mcEId].pt();
606					McProbe_rescEta_[iComb_] = McElecsResc[mcEId].eta();
607					McProbe_rescPhi_[iComb_] = McElecsResc[mcEId].phi();
608					probe_d_MCE_SCE_[iComb_] = McElecs[mcEId].energy() - it->second->superCluster()->rawEnergy();
609					McElecProbe_dPhi_[iComb_] = reco::deltaPhi(McElecs[mcEId].phi(), McElecs[(mcEId+1)%2].phi());
610					McElecProbe_dEta_[iComb_] = fabs(McElecs[mcEId].eta() - McElecs[(mcEId+1)%2].eta());
611					McElecProbe_dR_[iComb_] = reco::deltaR(McElecs[mcEId], McElecs[(mcEId+1)%2]);
612				} 
613			}
614
615			// Uncorrected supercluster V1
616			reco::SuperCluster scV1 = *(it->second->superCluster());
617			math::XYZTLorentzVector	probe_scV1_detVec = DetectorVector(scV1);
618			probe_sc_pt_[iComb_] = probe_scV1_detVec.pt(); 
619			probe_sc_eta_[iComb_] = scV1.eta();
620			probe_sc_phi_[iComb_] = scV1.phi();
621			probe_sc_nClus_[iComb_] = scV1.clustersSize();
622			probe_sc_E_[iComb_] = scV1.energy();
623			probe_sc_rawE_[iComb_] = scV1.rawEnergy();
624
625			ersatzMEt = ersatzFabrik(it->first, scV1, calomet, 1);
626			ErsatzV1CaloMEt_[iComb_] = ersatzMEt.pt();
627			ErsatzV1CaloMEtPhi_[iComb_] = ersatzMEt.phi();
628			//ersatzMEt = ersatzFabrik(it->first, it->second, t1met);
629			//ErsatzV1T1MEt_[iComb_] = ersatzMEt.pt();
630			//ErsatzV1T1MEtPhi_[iComb_] = ersatzMEt.phi();
631			ersatzMEt = ersatzFabrik(it->first, it->second, pfmet);
632			ErsatzV1PfMEt_[iComb_] = ersatzMEt.pt();
633			ErsatzV1PfMEtPhi_[iComb_] = ersatzMEt.phi();
634			ersatzMEt = ersatzFabrik(it->first, it->second, tcmet);
635			ErsatzV1TcMEt_[iComb_] = ersatzMEt.pt();
636			ErsatzV1TcMEtPhi_[iComb_] = ersatzMEt.phi();
637
638			// fEta corrected supercluster V2
639			reco::SuperCluster scV2;
640			if(fabs(probe_sc_eta_[iComb_]) < 1.479) 
641			{
642				scV2 = fEtaScCorr(scV1);
643			}else{
644				scV2 = scV1;
645			}
646			probe_scV2_E_[iComb_] = scV2.energy();
647			ersatzMEt = ersatzFabrik(it->first, scV2, calomet, 2);
648			ErsatzV2CaloMEt_[iComb_] = ersatzMEt.pt();
649			ErsatzV2CaloMEtPhi_[iComb_] = ersatzMEt.phi();
650		
651			// fBrem corrected supercluster V3
652			reco::SuperCluster scV3;	
653			if(fabs(probe_sc_eta_[iComb_]) < 1.479) 
654			{
655				scV3 = fBremScCorr(scV1, hyb_fCorrPSet_);
656			}else{
657				scV3 = fBremScCorr(scV1, m5x5_fCorrPSet_);
658			}	
659			probe_scV3_E_[iComb_] = scV3.energy();
660			ersatzMEt = ersatzFabrik(it->first, scV3, calomet, 3);
661			ErsatzV3CaloMEt_[iComb_] = ersatzMEt.pt();
662			ErsatzV3CaloMEtPhi_[iComb_] = ersatzMEt.phi();
663			
664			// Fully corrected supercluster V4
665			reco::SuperCluster scV4;
666			if(fabs(probe_sc_eta_[iComb_]) < 1.479) 
667			{
668				scV4 = fBremScCorr(scV1, hyb_fCorrPSet_);
669			}else{
670				scV4 = fBremScCorr(scV1, m5x5_fCorrPSet_);
671			}	
672			probe_scV4_E_[iComb_] = scV4.energy();
673			ersatzMEt = ersatzFabrik(it->first, scV4, calomet, 4);
674			ErsatzV4CaloMEt_[iComb_] = ersatzMEt.pt();
675			ErsatzV4CaloMEtPhi_[iComb_] = ersatzMEt.phi();
676
677			++iComb_;
678		}
679		t_->Fill();
680	}
681	}
682}
683
684std::map<reco::GsfElectronRef, reco::GsfElectronRef> ErsatzMEt::probeFinder(const std::vector<reco::GsfElectronRef>& tags,
685							const edm::Handle<reco::GsfElectronCollection> pElectrons)
686{
687	const reco::GsfElectronCollection *probeCands = pElectrons.product();
688	std::map<reco::GsfElectronRef, reco::GsfElectronRef> TagProbes;
689	for(std::vector<reco::GsfElectronRef>::const_iterator tagelec = tags.begin(); tagelec != tags.end(); ++tagelec)
690	{
691		reco::GsfElectronRef tag = *tagelec;
692		std::pair<reco::GsfElectronRef, reco::GsfElectronRef> TagProbePair;
693		int nProbesPerTag = 0;
694		int index = 0;
695		for(reco::GsfElectronCollection::const_iterator probeelec = probeCands->begin(); probeelec != probeCands->end(); ++probeelec)
696		{
697			reco::GsfElectronRef probe(pElectrons, index);
698			double probeScEta = probe->superCluster()->eta();
699			if(probe->superCluster() != tag->superCluster() && fabs(probeScEta) < 2.5)
700			{
701				if(fabs(probeScEta) < 1.4442 || fabs(probeScEta) > 1.560)
702				{
703					double invmass = ROOT::Math::VectorUtil::InvariantMass(tag->p4(), probe->p4());
704					if(mTPmin_ <= invmass && invmass <= mTPmax_)
705					{
706						TagProbePair = std::make_pair(tag, probe);
707						++nProbesPerTag;
708					}		
709				}
710			}
711			++index;
712		}
713		//nGsfElectrons_ = index;
714		if(nProbesPerTag == 1) TagProbes.insert(TagProbePair);
715	}	
716	return TagProbes;
717}
718
719reco::MET ErsatzMEt::ersatzFabrik(const reco::GsfElectronRef& elec, 
720					const reco::SuperCluster& sc,
721					const reco::MET& met,
722					const int corr)
723{
724	const math::XYZPoint ZVertex(elec->TrackPositionAtVtx().X(), elec->TrackPositionAtVtx().Y(),elec->TrackPositionAtVtx().Z());
725
726	math::XYZTLorentzVector nu, boost_nu, ele, boost_ele;
727	reco::SuperCluster elecSc = *(elec->superCluster());
728	nu = PhysicsVectorRaw(met.vertex(), sc);
729	boost_nu = PhysicsVectorRaw(ZVertex, sc);
730	ele = PhysicsVectorRaw(met.vertex(), elecSc);
731	boost_ele = ele;
732
733	//Should use reco vertex for best Z->ee measurement. 
734        edm::LogDebug_("ersatzFabrikV1", "", 569)<<"elec  = ("<< elec->p4().Px() << ", "<< elec->p4().Py()<< ", "<< elec->p4().Pz() << ", "<< elec->p4().E()<<")";
735	math::XYZTLorentzVector Zboson = boost_nu + elec->p4();
736        edm::LogDebug_("ersatzFabrikV1", "", 569)<<"Z pt = "<< Zboson.Pt() << "Z boson mass = " << Zboson.M();
737        edm::LogDebug_("ersatzFabrikV1","", 570)<<"Z boson in lab frame = ("<<Zboson.Px()<<", "<<Zboson.Py()<<", "
738                                        <<Zboson.Pz()<<", "<<Zboson.E()<<")";
739        math::XYZTLorentzVector RescZboson(Zboson.Px(), Zboson.Py(), Zboson.Pz(), sqrt(Zboson.P2()+(mW_*mW_*Zboson.M2())/(mZ_*mZ_)));
740        edm::LogDebug_("ersatzFabrikV1","", 573)<<"W boson in lab frame = ("<<RescZboson.Px()<<", "<<RescZboson.Py()<<", "
741                                        <<RescZboson.Pz()<<", "<<RescZboson.E()<<")";
742        ROOT::Math::Boost BoostToZRestFrame(Zboson.BoostToCM());
743        edm::LogDebug_("ersatzFabrikV1","", 576)<<"Electron in lab frame = ("<< boost_ele.Px()<<", "<< boost_ele.Py()<<", "
744                                        << boost_ele.Pz()<<", "<< boost_ele.E()<<")";
745        edm::LogDebug_("ersatzFabrikV1","", 578)<<"Ersatz Neutrino in lab frame = ("<< boost_nu.Px()<<", "<< boost_nu.Py()<<", "
746                                        << boost_nu.Pz()<<", "<< boost_nu.E()<<")";
747        boost_ele = BoostToZRestFrame(boost_ele);
748        boost_nu = BoostToZRestFrame(boost_nu);
749        edm::LogDebug_("ersatzFabrikV1","", 582)<<"Electron in Z rest frame = ("<<boost_ele.Px()<<", "<<boost_ele.Py()<<", "
750                                        <<boost_ele.Pz()<<", "<<boost_ele.E()<<")";
751        edm::LogDebug_("ersatzFabrikV1","", 584)<<"Ersatz Neutrino in Z rest frame = ("<<boost_nu.Px()<<", "<<boost_nu.Py()<<", "
752                                        <<boost_nu.Pz()<<", "<<boost_nu.E()<<")";
753        boost_ele *= mW_/mZ_;
754        boost_nu *= mW_/mZ_;
755
756        double E_W = RescZboson.E();
757        ROOT::Math::Boost BackToLab(RescZboson.Px()/E_W, RescZboson.Py()/E_W, RescZboson.Pz()/E_W);
758	math::XYZTLorentzVector metVec(-99999., -99999., -99., -99999.);
759        boost_ele = BackToLab(boost_ele);
760
761       	boost_nu = BackToLab(boost_nu);
762        math::XYZTLorentzVector sum = boost_nu+boost_ele;
763       	edm::LogDebug_("ersatzFabrikV1","", 597)<<"Electron back in lab frame = ("<<boost_ele.Px()<<", "<<boost_ele.Py()<<", "
764                                        <<boost_ele.Pz()<<", "<<boost_ele.E()<<")";
765        edm::LogDebug_("ersatzFabrikV1","", 599)<<"Ersatz Neutrino back in lab frame = ("<<boost_nu.Px()<<", "<<boost_nu.Py()<<", "
766                                       <<boost_nu.Pz()<<", "<<boost_nu.E()<<")";
767       	edm::LogDebug_("ersatzFabrikV1","", 601)<<"boost_ele + boost_nu = ("<<sum.Px()<<", "<<sum.Py()<<", "
768                                        <<sum.Pz()<<", "<<sum.E()<<")";
769
770	nu.SetXYZT(nu.X(), nu.Y(), 0., nu.T());
771	ele.SetXYZT(ele.X(), ele.Y(), 0., ele.T());
772	boost_ele.SetXYZT(boost_ele.X(), boost_ele.Y(), 0., boost_ele.T());
773	metVec = met.p4() + nu + ele - boost_ele;
774	
775	reco::MET ersatzMEt(metVec, met.vertex()); 	
776	if (corr == 1)	
777	{
778		//Z_caloV1_m_[iComb_] = Zboson.M();
779		//Z_caloV1_pt_[iComb_] = Zboson.Pt();
780		//Z_caloV1_y_[iComb_] = Zboson.Y();
781		//Z_caloV1_eta_[iComb_] = Zboson.Eta();
782		//Z_caloV1_phi_[iComb_] = Zboson.Phi();
783		//Z_caloV1_rescM_[iComb_] = RescZboson.M();
784		//Z_caloV1_rescPt_[iComb_] = RescZboson.Pt();
785		//Z_caloV1_rescY_[iComb_] = RescZboson.Y();
786		//Z_caloV1_rescEta_[iComb_] = RescZboson.Eta();
787		//Z_caloV1_rescPhi_[iComb_] = RescZboson.Phi();
788		//Z_caloV1_probe_dPhi_[iComb_] = reco::deltaPhi(Zboson.Phi(), elec->phi());
789		//tag_caloV1_rescPt_[iComb_]  = boost_ele.Pt();
790		//tag_caloV1_rescEta_[iComb_]  = boost_ele.Eta();
791		//tag_caloV1_rescPhi_[iComb_]  = boost_ele.Phi();
792		//probe_caloV1_rescPt_[iComb_]  = boost_nu.Pt();
793		//probe_caloV1_rescEta_[iComb_]  = boost_nu.Eta();
794		//probe_caloV1_rescPhi_[iComb_]  = boost_nu.Phi();
795		ErsatzV1_Mesc_[iComb_]  = ROOT::Math::VectorUtil::InvariantMass(elec->p4(), boost_nu);
796		ErsatzV1_rescMesc_[iComb_] = ROOT::Math::VectorUtil::InvariantMass(ele, nu);
797		ErsatzV1CaloMt_[iComb_] = sqrt(2.*boost_ele.Pt()*ersatzMEt.pt()*
798						(1-cos(reco::deltaPhi(boost_ele.Phi(), ersatzMEt.phi()))));
799	}
800	if (corr == 2) 	
801	{
802		//Z_caloV2_m_[iComb_] = Zboson.M();
803		//Z_caloV2_pt_[iComb_] = Zboson.Pt();
804		//Z_caloV2_y_[iComb_] = Zboson.Y();
805		//Z_caloV2_eta_[iComb_] = Zboson.Eta();
806		//Z_caloV2_phi_[iComb_] = Zboson.Phi();
807		//Z_caloV2_rescM_[iComb_] = RescZboson.M();
808		//Z_caloV2_rescPt_[iComb_] = RescZboson.Pt();
809		//Z_caloV2_rescY_[iComb_] = RescZboson.Y();
810		//Z_caloV2_rescEta_[iComb_] = RescZboson.Eta();
811		//Z_caloV2_rescPhi_[iComb_] = RescZboson.Phi();
812		//Z_caloV2_probe_dPhi_[iComb_] = reco::deltaPhi(Zboson.Phi(), boost_elec->phi());
813		//tag_caloV2_rescPt_[iComb_]  = boost_ele.Pt();
814		//tag_caloV2_rescEta_[iComb_]  = boost_ele.Eta();
815		//tag_caloV2_rescPhi_[iComb_]  = boost_ele.Phi();
816		//probe_caloV2_rescPt_[iComb_]  = boost_nu.Pt();
817		//probe_caloV2_rescEta_[iComb_]  = boost_nu.Eta();
818		//probe_caloV2_rescPhi_[iComb_]  = boost_nu.Phi();
819		ErsatzV2_Mesc_[iComb_]  = ROOT::Math::VectorUtil::InvariantMass(elec->p4(), boost_nu);
820		ErsatzV2_rescMesc_[iComb_] = ROOT::Math::VectorUtil::InvariantMass(ele, nu);
821		ErsatzV2CaloMt_[iComb_] = sqrt(2.*boost_ele.Pt()*ersatzMEt.pt()*
822						(1-cos(reco::deltaPhi(boost_ele.Phi(), ersatzMEt.phi()))));
823	}
824	if (corr == 3) 	
825	{
826		//Z_caloV3_m_[iComb_] = Zboson.M();
827		//Z_caloV3_pt_[iComb_] = Zboson.Pt();
828		//Z_caloV3_y_[iComb_] = Zboson.Y();
829		//Z_caloV3_eta_[iComb_] = Zboson.Eta();
830		//Z_caloV3_phi_[iComb_] = Zboson.Phi();
831		//Z_caloV3_rescM_[iComb_] = RescZboson.M();
832		//Z_caloV3_rescPt_[iComb_] = RescZboson.Pt();
833		//Z_caloV3_rescY_[iComb_] = RescZboson.Y();
834		//Z_caloV3_rescEta_[iComb_] = RescZboson.Eta();
835		//Z_caloV3_rescPhi_[iComb_] = RescZboson.Phi();
836		//Z_caloV3_probe_dPhi_[iComb_] = reco::deltaPhi(Zboson.Phi(), boost_elec->phi());
837		//tag_caloV3_rescPt_[iComb_]  = boost_ele.Pt();
838		//tag_caloV3_rescEta_[iComb_]  = boost_ele.Eta();
839		//tag_caloV3_rescPhi_[iComb_]  = boost_ele.Phi();
840		//probe_caloV3_rescPt_[iComb_]  = boost_nu.Pt();
841		//probe_caloV3_rescEta_[iComb_]  = boost_nu.Eta();
842		//probe_caloV3_rescPhi_[iComb_]  = boost_nu.Phi();
843		ErsatzV3_Mesc_[iComb_]  = ROOT::Math::VectorUtil::InvariantMass(elec->p4(), boost_nu);
844		ErsatzV3_rescMesc_[iComb_] = ROOT::Math::VectorUtil::InvariantMass(ele, nu);
845		ErsatzV3CaloMt_[iComb_] = sqrt(2.*boost_ele.Pt()*ersatzMEt.pt()*
846						(1-cos(reco::deltaPhi(boost_ele.Phi(), ersatzMEt.phi()))));
847	}
848	if (corr == 4) 	
849	{
850		//Z_caloV4_m_[iComb_] = Zboson.M();
851		//Z_caloV4_pt_[iComb_] = Zboson.Pt();
852		//Z_caloV4_y_[iComb_] = Zboson.Y();
853		//Z_caloV4_eta_[iComb_] = Zboson.Eta();
854		//Z_caloV4_phi_[iComb_] = Zboson.Phi();
855		//Z_caloV4_rescM_[iComb_] = RescZboson.M();
856		//Z_caloV4_rescPt_[iComb_] = RescZboson.Pt();
857		//Z_caloV4_rescY_[iComb_] = RescZboson.Y();
858		//Z_caloV4_rescEta_[iComb_] = RescZboson.Eta();
859		//Z_caloV4_rescPhi_[iComb_] = RescZboson.Phi();
860		//Z_caloV4_probe_dPhi_[iComb_] = reco::deltaPhi(Zboson.Phi(), boost_elec->phi());
861		//tag_caloV4_rescPt_[iComb_]  = boost_ele.Pt();
862		//tag_caloV4_rescEta_[iComb_]  = boost_ele.Eta();
863		//tag_caloV4_rescPhi_[iComb_]  = boost_ele.Phi();
864		//probe_caloV4_rescPt_[iComb_]  = boost_nu.Pt();
865		//probe_caloV4_rescEta_[iComb_]  = boost_nu.Eta();
866		//probe_caloV4_rescPhi_[iComb_]  = boost_nu.Phi();
867		ErsatzV4_Mesc_[iComb_]  = ROOT::Math::VectorUtil::InvariantMass(elec->p4(), boost_nu);
868		ErsatzV4_rescMesc_[iComb_] = ROOT::Math::VectorUtil::InvariantMass(ele, nu);
869		ErsatzV4CaloMt_[iComb_] = sqrt(2.*boost_ele.Pt()*ersatzMEt.pt()*
870						(1-cos(reco::deltaPhi(boost_ele.Phi(), ersatzMEt.phi()))));
871	}
872	return ersatzMEt;
873}
874
875reco::MET ErsatzMEt::ersatzFabrik(const reco::GsfElectronRef& tag, 
876					const reco::GsfElectronRef& probe,
877					const reco::MET& met)
878{
879	math::XYZTLorentzVector elec, nu, boost_elec, boost_nu;
880	boost_elec = tag->p4();
881        edm::LogDebug_("ersatzFabrikV1", "", 858)<<"boost_elec  = ("<< boost_elec.Px() << ", "<< boost_elec.Py()<< ", "<< boost_elec.Pz() << ", "<< boost_elec.E()<<")";
882	boost_nu = probe->p4();
883        edm::LogDebug_("ersatzFabrikV1", "", 860)<<"boost_nu  = ("<< boost_nu.Px() << ", "<< boost_nu.Py()<< ", "<< boost_nu.Pz() << ", "<< boost_nu.E()<<")";
884	math::XYZTLorentzVector Zboson = boost_elec + boost_nu;
885        edm::LogDebug_("ersatzFabrikV1", "", 862)<<"Zboson  = ("<< Zboson.Px() << ", "<< Zboson.Py()<< ", "<< Zboson.Pz() << ", "<< Zboson.E()<<")";
886        math::XYZTLorentzVector RescZboson(Zboson.Px(), Zboson.Py(), Zboson.Pz(), sqrt(Zboson.P2()+(mW_*mW_*Zboson.M2())/(mZ_*mZ_)));
887        edm::LogDebug_("ersatzFabrikV1", "", 864)<<"RescZboson  = ("<< RescZboson.Px() << ", "<< RescZboson.Py()<< ", "<< RescZboson.Pz() << ", "<< RescZboson.E()<<")";
888        ROOT::Math::Boost BoostToZRestFrame(Zboson.BoostToCM());
889        elec = BoostToZRestFrame(boost_elec);
890        edm::LogDebug_("ersatzFabrikV1", "", 867)<<"boost_elec (in Z rest frame) = ("<< elec.Px() << ", "<< elec.Py()<< ", "<< elec.Pz() << ", "<< elec.E()<<")";
891        nu = BoostToZRestFrame(boost_nu);
892        edm::LogDebug_("ersatzFabrikV1", "", 869)<<"boost_nu (in Z rest frame) = ("<< nu.Px() << ", "<< nu.Py()<< ", "<< nu.Pz() << ", "<< nu.E()<<")";
893        elec *= mW_/mZ_;
894        edm::LogDebug_("ersatzFabrikV1", "", 871)<<"elec (in Z rest frame) = ("<< elec.Px() << ", "<< elec.Py()<< ", "<< elec.Pz() << ", "<< elec.E()<<")";
895        nu *= mW_/mZ_;
896        edm::LogDebug_("ersatzFabrikV1", "", 873)<<"nu (in Z rest frame) = ("<< nu.Px() << ", "<< nu.Py()<< ", "<< nu.Pz() << ", "<< nu.E()<<")";
897        ROOT::Math::Boost BoostBackToLab(RescZboson.Px()/RescZboson.E(), RescZboson.Py()/RescZboson.E(), RescZboson.Pz()/RescZboson.E());
898	math::XYZTLorentzVector metVec(-99999., -99999., -99., -99999.);
899        elec = BoostBackToLab(elec);
900        edm::LogDebug_("ersatzFabrikV1", "", 877)<<"elec = ("<< elec.Px() << ", "<< elec.Py()<< ", "<< elec.Pz() << ", "<< elec.E()<<")";
901       	nu = BoostBackToLab(nu);
902        edm::LogDebug_("ersatzFabrikV1", "", 879)<<"nu = ("<< nu.Px() << ", "<< nu.Py()<< ", "<< nu.Pz() << ", "<< nu.E()<<")";
903	Z_m_[iComb_] = Zboson.M();
904	Z_pt_[iComb_] = Zboson.Pt();
905	Z_y_[iComb_] = Zboson.Y();
906	Z_eta_[iComb_] = Zboson.Eta();
907	Z_phi_[iComb_] = Zboson.Phi();
908	Z_rescM_[iComb_] = RescZboson.M();
909	Z_rescPt_[iComb_] = RescZboson.Pt();
910	Z_rescY_[iComb_] = RescZboson.Y();
911	Z_rescEta_[iComb_] = RescZboson.Eta();
912	Z_rescPhi_[iComb_] = RescZboson.Phi();
913	Z_probe_dPhi_[iComb_] = reco::deltaPhi(Zboson.Phi(), boost_elec.phi());
914	tag_rescPt_[iComb_]  = elec.Pt();
915	tag_rescEta_[iComb_]  = elec.Eta();
916	tag_rescPhi_[iComb_]  = elec.Phi();
917	probe_rescPt_[iComb_]  = nu.Pt();
918	probe_rescEta_[iComb_]  = nu.Eta();
919	probe_rescPhi_[iComb_]  = nu.Phi();
920	elec.SetXYZT(elec.X(), elec.Y(), 0., elec.T());
921	nu.SetXYZT(nu.X(), nu.Y(), 0., nu.T());
922	boost_elec.SetXYZT(boost_elec.X(), boost_elec.Y(), 0., boost_elec.T());
923	metVec = met.p4() + nu + elec - boost_elec;
924	reco::MET ersatzMEt(metVec, met.vertex()); 	
925	return ersatzMEt;
926}
927
928bool ErsatzMEt::isInBarrel(double eta)
929{
930	return (fabs(eta) < BarrelEtaMax_);
931}
932
933bool ErsatzMEt::isInEndCap(double eta)
934{
935	return (fabs(eta) < EndCapEtaMax_ && fabs(eta) > EndCapEtaMin_);
936}
937
938bool ErsatzMEt::isInFiducial(double eta)
939{
940	return isInBarrel(eta) || isInEndCap(eta);
941}
942
943// ------------ method called once each job just after ending the event loop  ------------
944void ErsatzMEt::endJob() {
945}
946//define this as a plug-in
947DEFINE_FWK_MODULE(ErsatzMEt);
948