/ElectroWeakAnalysis/ZEE/src/ErsatzMEt.cc

https://github.com/aivanov-cern/cmssw · C++ · 948 lines · 777 code · 68 blank · 103 comment · 61 complexity · 2313a34fb4c823d43c26c4f86f2bc9f0 MD5 · raw file

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