CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_0/src/ElectroWeakAnalysis/ZEE/src/ErsatzMEt.cc

Go to the documentation of this file.
00001 #include "ElectroWeakAnalysis/ZEE/interface/ErsatzMEt.h"
00002 #include "FWCore/Common/interface/TriggerNames.h"
00003 
00004 ErsatzMEt::ErsatzMEt(const edm::ParameterSet& ps)
00005 {
00006         MCTruthCollection_  = ps.getParameter<edm::InputTag>("MCTruthCollection");
00007         ElectronCollection_ = ps.getParameter<edm::InputTag>("ElectronCollection");
00008         HybridScCollection_ = ps.getParameter<edm::InputTag>("HybridScCollection");
00009         M5x5ScCollection_ = ps.getParameter<edm::InputTag>("M5x5ScCollection");
00010         GenMEtCollection_  = ps.getParameter<edm::InputTag>("GenMEtCollection");
00011         CaloMEtCollection_ = ps.getParameter<edm::InputTag>("CaloMEtCollection");
00012         //T1MEtCollection_ = ps.getParameter<edm::InputTag>("T1MEtCollection");
00013         PfMEtCollection_ = ps.getParameter<edm::InputTag>("PfMEtCollection");
00014         TcMEtCollection_ = ps.getParameter<edm::InputTag>("TcMEtCollection");
00015         TriggerEvent_ = ps.getParameter<edm::InputTag>("TriggerEvent");
00016         TriggerPath_ = ps.getParameter<edm::InputTag>("TriggerPath");
00017         TriggerResults_ = ps.getParameter<edm::InputTag>("TriggerResults");
00018         TriggerName_ = ps.getParameter<std::string>("TriggerName");
00019         HLTPathCheck_ = ps.getParameter<bool>("HLTPathCheck");
00020         
00021         Zevent_ = ps.getParameter<bool>("Zevent");
00022         mW_ = ps.getParameter<double>("mW");
00023         mZ_ = ps.getParameter<double>("mZ");
00024         mTPmin_ = ps.getParameter<double>("mTPmin");
00025         mTPmax_ = ps.getParameter<double>("mTPmax");
00026         BarrelEtaMax_ = ps.getParameter<double>("BarrelEtaMax");
00027         EndCapEtaMin_ = ps.getParameter<double>("EndCapEtaMin");
00028         EndCapEtaMax_ = ps.getParameter<double>("EndCapEtaMax");
00029         
00030         hyb_fCorrPSet_ = ps.getParameter<edm::ParameterSet>("hyb_fCorrPSet");
00031         m5x5_fCorrPSet_ = ps.getParameter<edm::ParameterSet>("m5x5_fCorrPSet");
00032 
00033         double CElecPtMin = ps.getParameter<double>("CElecPtMin");
00034         double CEB_siEiE = ps.getParameter<double>("CEB_sigmaIEtaIEta"); 
00035         double CEB_dPhiIn = ps.getParameter<double>("CEB_deltaPhiIn"); 
00036         double CEB_dEtaIn = ps.getParameter<double>("CEB_deltaEtaIn"); 
00037         double CEB_EcalIso = ps.getParameter<double>("CEB_EcalIso");
00038         double CEB_HcalIso = ps.getParameter<double>("CEB_HcalIso");
00039         double CEB_TrckIso = ps.getParameter<double>("CEB_TrckIso");
00040         double CEE_siEiE = ps.getParameter<double>("CEE_sigmaIEtaIEta"); 
00041         double CEE_dPhiIn = ps.getParameter<double>("CEE_deltaPhiIn"); 
00042         double CEE_dEtaIn = ps.getParameter<double>("CEE_deltaEtaIn"); 
00043         double CEE_EcalIso = ps.getParameter<double>("CEE_EcalIso");
00044         double CEE_HcalIso = ps.getParameter<double>("CEE_HcalIso");
00045         double CEE_TrckIso = ps.getParameter<double>("CEE_TrckIso");
00046 
00047         CutVector_.resize(13);
00048         CutVector_[EtCut_] = CElecPtMin;
00049         CutVector_[EB_sIhIh_] = CEB_siEiE;
00050         CutVector_[EB_dPhiIn_] = CEB_dPhiIn;
00051         CutVector_[EB_dEtaIn_] = CEB_dEtaIn;
00052         CutVector_[EB_TrckIso_] = CEB_TrckIso;
00053         CutVector_[EB_EcalIso_] = CEB_EcalIso;
00054         CutVector_[EB_HcalIso_] = CEB_HcalIso;
00055         CutVector_[EE_sIhIh_] = CEE_siEiE;
00056         CutVector_[EE_dPhiIn_] = CEE_dPhiIn;
00057         CutVector_[EE_dEtaIn_] = CEE_dEtaIn;
00058         CutVector_[EE_TrckIso_] = CEE_TrckIso;
00059         CutVector_[EE_EcalIso_] = CEE_EcalIso;
00060         CutVector_[EE_HcalIso_] = CEE_HcalIso;
00061         
00062         for(std::vector<double>::const_iterator it = CutVector_.begin(); it != CutVector_.end(); ++it)
00063         {
00064                 edm::LogDebug_("","",101)<<"CutVector_ = "<< *it;
00065         }
00066 }
00067 
00068 ErsatzMEt::~ErsatzMEt()
00069 {
00070 }
00071 
00072 // ------------ method called once each job just before starting event loop  ------------
00073 void ErsatzMEt::beginJob()
00074 {
00075         edm::Service<TFileService> fs;
00076 
00077         t_ = fs->make<TTree>("ErsatzMEt", "Data on ErsatzMEt");
00078         
00079         edm::LogDebug_("","", 75)<<"Creating Ersatz MEt branches.";
00080         
00081         t_->Branch("nTags", &nTags_, "nTags/I");
00082         t_->Branch("nProbes", &nProbes_, "nProbes/I");
00083         
00084         t_->Branch("ErsatzV1CaloMEt", ErsatzV1CaloMEt_, "ErsatzV1CaloMEt[4]/D");
00085         t_->Branch("ErsatzV1CaloMt", ErsatzV1CaloMt_, "ErsatzV1CaloMt[4]/D");
00086         t_->Branch("ErsatzV1CaloMEtPhi", ErsatzV1CaloMEtPhi_, "ErsatzV1CaloMEtPhi[4]/D");
00087         t_->Branch("ErsatzV2CaloMEt", ErsatzV2CaloMEt_, "ErsatzV2CaloMEt[4]/D");
00088         t_->Branch("ErsatzV2CaloMEtPhi", ErsatzV2CaloMEtPhi_, "ErsatzV2CaloMEtPhi[4]/D");
00089         t_->Branch("ErsatzV2CaloMt", ErsatzV2CaloMt_, "ErsatzV2CaloMt[4]/D");
00090         t_->Branch("ErsatzV3CaloMEt", ErsatzV3CaloMEt_, "ErsatzV3CaloMEt[4]/D");
00091         t_->Branch("ErsatzV3CaloMEtPhi", ErsatzV3CaloMEtPhi_, "ErsatzV3CaloMEtPhi[4]/D");
00092         t_->Branch("ErsatzV3CaloMt", ErsatzV3CaloMt_, "ErsatzV3CaloMt[4]/D");
00093         t_->Branch("ErsatzV4CaloMEt", ErsatzV4CaloMEt_, "ErsatzV4CaloMEt[4]/D");
00094         t_->Branch("ErsatzV4CaloMEtPhi", ErsatzV4CaloMEtPhi_, "ErsatzV4CaloMEtPhi[4]/D");
00095         t_->Branch("ErsatzV4CaloMt", ErsatzV4CaloMt_, "ErsatzV4CaloMt[4]/D");
00096         t_->Branch("ErsatzV1T1MEt", ErsatzV1T1MEt_, "ErsatzV1T1MEt[4]/D");
00097         t_->Branch("ErsatzV1T1Mt", ErsatzV1T1Mt_, "ErsatzV1T1Mt[4]/D");
00098         t_->Branch("ErsatzV1T1MEtPhi", ErsatzV1T1MEtPhi_, "ErsatzV1T1MEtPhi[4]/D");
00099         t_->Branch("ErsatzV1PfMEt", ErsatzV1PfMEt_, "ErsatzV1PfMEt[4]/D");
00100         t_->Branch("ErsatzV1PfMt", ErsatzV1PfMt_, "ErsatzV1PfMt[4]/D");
00101         t_->Branch("ErsatzV1PfMEtPhi", ErsatzV1TcMEtPhi_, "ErsatzV1PfMEtPhi[4]/D");
00102         t_->Branch("ErsatzV1TcMEt", ErsatzV1TcMEt_, "ErsatzV1TcMEt[4]/D");
00103         t_->Branch("ErsatzV1TcMt", ErsatzV1TcMt_, "ErsatzV1TcMt[4]/D");
00104         t_->Branch("ErsatzV1TcMEtPhi", ErsatzV1TcMEtPhi_, "ErsatzV1TcMEtPhi[4]/D");
00105 
00106         t_->Branch("CaloMEt", &CaloMEt_, "CaloMEt/D");
00107         t_->Branch("CaloMEtphi", &CaloMEtphi_, "CaloMEtphi/D");
00108         t_->Branch("T1MEt", &T1MEt_, "T1MEt/D");
00109         t_->Branch("T1MEtphi", &T1MEtphi_, "T1MEtphi/D");
00110         t_->Branch("PfMEt", &PfMEt_, "PfMEt/D");
00111         t_->Branch("PfMEtphi", &PfMEtphi_, "PfMEtphi/D");
00112         t_->Branch("TcMEt", &TcMEt_, "TcMEt/D");
00113         t_->Branch("TcMEtphi", &TcMEtphi_, "TcMEtphi/D");
00114 
00115         edm::LogDebug_("","", 91)<<"Creating electron branches.";
00116         t_->Branch("tag_q", tag_q_,"tag_q[4]/I");
00117         t_->Branch("tag_pt", tag_pt_,"tag_pt[4]/D");
00118         t_->Branch("tag_eta", tag_eta_,"tag_eta[4]/D");
00119         t_->Branch("tag_phi", tag_phi_,"tag_phi[4]/D");
00120         t_->Branch("tag_sIhIh", tag_sIhIh_, "tag_sIhIh[4]/D");
00121         t_->Branch("tag_dPhiIn", tag_dPhiIn_, "tag_dPhiIn[4]/D");
00122         t_->Branch("tag_dEtaIn", tag_dEtaIn_, "tag_dEtaIn[4]/D");
00123         t_->Branch("tag_trckIso", tag_trckIso_,"tag_trckIso[4]/D");
00124         t_->Branch("tag_ecalIso", tag_ecalIso_,"tag_ecalIso[4]/D");
00125         t_->Branch("tag_hcalIso", tag_hcalIso_,"tag_hcalIso[4]/D");
00126         t_->Branch("tag_e2x5Max", tag_e2x5Max_,"tag_e2x5Max[4]/D");
00127         t_->Branch("tag_e1x5Max", tag_e1x5Max_,"tag_e1x5Max[4]/D");
00128         t_->Branch("tag_e5x5", tag_e5x5_,"tag_e5x5[4]/D");
00129         t_->Branch("tag_hoe", tag_hoe_,"tag_hoe[4]/D");
00130         t_->Branch("tag_eop", tag_eop_,"tag_eop[4]/D");
00131         t_->Branch("tag_pin", tag_pin_,"tag_pin[4]/D");
00132         t_->Branch("tag_pout", tag_pout_,"tag_pout[4]/D");
00133         t_->Branch("tag_rescPt", tag_rescPt_, "tag_rescPt[4]/D");
00134         t_->Branch("tag_rescEta", tag_rescEta_, "tag_rescEta[4]/D");
00135         t_->Branch("tag_rescPhi", tag_rescPhi_, "tag_rescPhi[4]/D");
00136         
00137         edm::LogDebug_("","", 103)<<"Creating ersatz neutrino branches.";
00138         t_->Branch("probe_q", probe_q_,"probe_q[4]/I");
00139         t_->Branch("probe_pt", probe_pt_,"probe_pt[4]/D");
00140         t_->Branch("probe_eta", probe_eta_,"probe_eta[4]/D");
00141         t_->Branch("probe_phi", probe_phi_,"probe_phi[4]/D");
00142         t_->Branch("probe_sIhIh", probe_sIhIh_, "probe_sIhIh[4]/D");
00143         t_->Branch("probe_dPhiIn", probe_dPhiIn_, "probe_dPhiIn[4]/D");
00144         t_->Branch("probe_dEtaIn", probe_dEtaIn_, "probe_dEtaIn[4]/D");
00145         t_->Branch("probe_trckIso", probe_trckIso_,"probe_trckIso[4]/D");
00146         t_->Branch("probe_ecalIso", probe_ecalIso_,"probe_ecalIso[4]/D");
00147         t_->Branch("probe_hcalIso", probe_hcalIso_,"probe_hcalIso[4]/D");
00148         t_->Branch("probe_e2x5Max", probe_e2x5Max_,"probe_e2x5Max[4]/D");
00149         t_->Branch("probe_e1x5Max", probe_e1x5Max_,"probe_e1x5Max[4]/D");
00150         t_->Branch("probe_e5x5", probe_e5x5_,"probe_e5x5[4]/D");
00151         t_->Branch("probe_hoe", probe_hoe_,"probe_hoe[4]/D");
00152         t_->Branch("probe_eop", probe_eop_,"probe_eop[4]/D");
00153         t_->Branch("probe_pin", probe_pin_,"probe_pin[4]/D");
00154         t_->Branch("probe_pout", probe_pout_,"probe_pout[4]/D");
00155         t_->Branch("probe_rescPt", probe_rescPt_, "probe_rescPt[4]/D");
00156         t_->Branch("probe_rescEta", probe_rescEta_, "probe_rescEta[4]/D");
00157         t_->Branch("probe_rescPhi", probe_rescPhi_, "probe_rescPhi[4]/D");
00158 
00159         t_->Branch("Z_m", Z_m_, "Z_m[4]/D");
00160         t_->Branch("Z_pt", Z_pt_, "Z_pt[4]/D");
00161         t_->Branch("Z_eta", Z_eta_, "Z_eta[4]/D");
00162         t_->Branch("Z_y", Z_y_, "Z_y[4]/D");
00163         t_->Branch("Z_phi", Z_phi_, "Z_phi[4]/D");
00164         t_->Branch("Z_rescM", Z_rescM_, "Z_rescM[4]/D");
00165         t_->Branch("Z_rescPt", Z_rescPt_, "Z_rescPt[4]/D");
00166         t_->Branch("Z_rescEta", Z_rescEta_, "Z_rescEta[4]/D");
00167         t_->Branch("Z_rescY", Z_rescY_, "Z_rescY[4]/D");
00168         t_->Branch("Z_rescPhi", Z_rescPhi_, "Z_rescPhi[4]/D");
00169         t_->Branch("Z_probe_dPhi",Z_probe_dPhi_,"Z_probe_dPhi[4]/D"); 
00170 
00171         t_->Branch("probe_sc_pt",probe_sc_pt_, "probe_sc_pt[4]/D");
00172         t_->Branch("probe_sc_eta",probe_sc_eta_, "probe_sc_eta[4]/D");
00173         t_->Branch("probe_sc_phi", probe_sc_phi_, "probe_sc_phi[4]/D");
00174         t_->Branch("probe_sc_E",probe_sc_E_, "probe_sc_E[4]/D");
00175         t_->Branch("probe_sc_rawE",probe_sc_rawE_, "probe_sc_rawE[4]/D");
00176         t_->Branch("probe_sc_nClus", probe_sc_nClus_, "probe_sc_nClus[4]/D");
00177         t_->Branch("probe_scV2_E",probe_scV2_E_, "probe_scV2_E[4]/D");
00178         t_->Branch("probe_scV3_E",probe_scV3_E_, "probe_scV3_E[4]/D");
00179         t_->Branch("probe_scV4_E",probe_scV4_E_, "probe_scV4_E[4]/D");
00180         t_->Branch("probe_d_MCE_SCE", probe_d_MCE_SCE_, "probe_d_MCE_SCE[4]/D");
00181 
00182         t_->Branch("ErsatzV1_Mesc", ErsatzV1_Mesc_, "ErsatzV1_Mesc[4]/D");
00183         t_->Branch("ErsatzV1_rescMesc", ErsatzV1_rescMesc_, "ErsatzV1_rescMesc[4]/D");
00184 
00185         t_->Branch("McElec_nFinal", &McElec_nFinal_, "McElec_nFinal/I");
00186 
00187         if(Zevent_){
00188                 t_->Branch("McZ_m", &McZ_m_, "McZ_m/D");
00189                 t_->Branch("McZ_rescM", &McZ_rescM_, "McZ_rescM/D");
00190                 t_->Branch("McZ_Pt", &McZ_pt_, "McZ_Pt/D");
00191                 t_->Branch("McZ_y", &McZ_y_, "McZ_y/D");
00192                 t_->Branch("McZ_Eta", &McZ_eta_, "McZ_Eta/D");
00193                 t_->Branch("McZ_Phi", &McZ_phi_, "McZ_Phi/D");
00194                 t_->Branch("McZ_rescPt", &McZ_rescPt_, "McZ_Pt/D");
00195                 t_->Branch("McZ_rescY", &McZ_rescY_, "McZ_rescY/D");
00196                 t_->Branch("McZ_rescEta", &McZ_rescEta_, "McZ_Eta/D");
00197                 t_->Branch("McZ_rescPhi", &McZ_rescPhi_, "McZ_Phi/D");
00198                 t_->Branch("McElec_nZmum", &McElec_nZmum_, "McElec_nZmum/I");
00199                 t_->Branch("McElec_eta", &McElec_eta_, "McElec_eta[4]/D");
00200                 t_->Branch("McElec_pt", &McElec_pt_, "McElec_pt[4]/D");
00201                 t_->Branch("McElec_phi", &McElec_phi_, "McElec_phi[4]/D");
00202                 t_->Branch("McElec_rescEta", &McElec_rescEta_, "McElec_rescEta[4]/D");
00203                 t_->Branch("McElec_rescPhi", &McElec_rescPhi_, "McElec_rescPhi[4]/D");
00204                 t_->Branch("McElec_rescPt", &McElec_rescPt_, "McElec_rescPt[4]/D");
00205                 t_->Branch("McProbe_eta", &McProbe_eta_, "McProbe_eta[4]/D");
00206                 t_->Branch("McProbe_pt", &McProbe_pt_, "McProbe_pt[4]/D");
00207                 t_->Branch("McProbe_phi", &McProbe_phi_, "McProbe_phi[4]/D");
00208                 t_->Branch("McProbe_rescEta", &McProbe_rescEta_, "McProbe_rescEta[4]/D");
00209                 t_->Branch("McProbe_rescPt", &McProbe_rescPt_, "McProbe_rescPt[4]/D");
00210                 t_->Branch("McProbe_rescPhi", &McProbe_rescPhi_, "McProbe_rescPhi[4]/D");
00211                 t_->Branch("McElecProbe_dPhi", &McElecProbe_dPhi_, "McElecProbe_dPhi/D");
00212                 t_->Branch("McElecProbe_dEta", &McElecProbe_dEta_, "McElecProbe_dEta/D");
00213                 t_->Branch("McElecProbe_dR", &McElecProbe_dR_, "McElecProbe_dR/D");
00214         }
00215 
00216 }
00217 
00218 void ErsatzMEt::analyze(const edm::Event& evt, const edm::EventSetup& es)
00219 {
00220 
00221        es.get<CaloGeometryRecord>().get(geoHandle_);
00222        es.get<CaloTopologyRecord>().get(pTopology_);
00223 
00224         edm::LogDebug_("","", 151)<<"Initialising variables.";
00225         nTags_ = -99; nProbes_ = -99;
00226         CaloMEt_ = -99.; CaloMEtphi_ = -99.;
00227         T1MEt_ = -99.; T1MEtphi_ = -99.;
00228         PfMEt_ = -99.; PfMEtphi_ = -99.;
00229         TcMEt_ = -99.; TcMEtphi_ = -99.;        
00230         if(Zevent_)
00231         {
00232                 McZ_m_ = -99.; McZ_pt_ = -99.; McZ_y_ = -99.; McZ_eta_ = -99.; McZ_phi_ = -99.;
00233                 McZ_rescM_ = -99.; McZ_rescPt_ = -99.; McZ_rescY_ = -99.; McZ_rescEta_ = -99.; McZ_rescPhi_ = -99.;
00234                 McElec_nZmum_ = -99; McElec_nFinal_ = -99;
00235         }
00236 
00237         for(int i = 0; i < nEntries_arr_; ++i)
00238         {
00239                 tag_q_[i] = -99;
00240                 tag_pt_[i] = -99.; tag_eta_[i] = -99.; tag_phi_[i] = -99.;
00241                 tag_rescPt_[i] = -99.; tag_rescEta_[i] = -99.; tag_rescPhi_[i] = -99.;
00242                 tag_trckIso_[i] = -99.; tag_ecalIso_[i] = -99.; tag_hcalIso_[i] = -99.;
00243                 tag_sIhIh_[i] = -99.; tag_dPhiIn_[i] = -99.; tag_dEtaIn_[i] = -99.;
00244                 tag_e5x5_[i] = -99.; tag_e2x5Max_[i] = -99.; tag_e1x5Max_[i] = -99.;
00245                 tag_hoe_[i] = -99.; tag_eop_[i] = -99.; tag_pin_[i] = -99.; tag_pout_[i] = -99.;
00246 
00247                 probe_q_[i] = -99;
00248                 probe_pt_[i] = -99.; probe_eta_[i] = -99.; probe_phi_[i] = -99.;
00249                 probe_rescPt_[i] = -99.; probe_rescEta_[i] = -99.; probe_rescPhi_[i] = -99.;
00250                 probe_trckIso_[i] = -99.; probe_ecalIso_[i] = -99.; probe_hcalIso_[i] = -99.;
00251                 probe_sIhIh_[i] = -99.; probe_dPhiIn_[i] = -99.; probe_dEtaIn_[i] = -99.;
00252                 probe_e5x5_[i] = -99.; probe_e2x5Max_[i] = -99.; probe_e1x5Max_[i] = -99.;
00253                 probe_hoe_[i] = -99.; probe_eop_[i] = -99.; probe_pin_[i] = -99.; probe_pout_[i] = -99.;
00254                 
00255                 Z_pt_[i] = -99.; Z_y_[i] = -99.; Z_eta_[i] = -99.; Z_phi_[i] = -99.; Z_m_[i] = -99.;
00256                 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.;
00257                 
00258                 ErsatzV1_Mesc_[i] = -99.; ErsatzV1_rescMesc_[i] = -99.;
00259                 ErsatzV2_Mesc_[i] = -99.; ErsatzV2_rescMesc_[i] = -99.;
00260                 ErsatzV3_Mesc_[i] = -99.; ErsatzV3_rescMesc_[i] = -99.;
00261                 ErsatzV4_Mesc_[i] = -99.; ErsatzV4_rescMesc_[i] = -99.;
00262                 ErsatzV1CaloMEt_[i] = -99.; ErsatzV1CaloMt_[i] = -99.; ErsatzV1CaloMEtPhi_[i] = -99.; 
00263                 ErsatzV2CaloMEt_[i] = -99.; ErsatzV2CaloMt_[i] = -99.; ErsatzV2CaloMEtPhi_[i] = -99.; 
00264                 ErsatzV3CaloMEt_[i] = -99.; ErsatzV3CaloMt_[i] = -99.; ErsatzV3CaloMEtPhi_[i] = -99.; 
00265                 ErsatzV4CaloMEt_[i] = -99.; ErsatzV4CaloMt_[i] = -99.; ErsatzV4CaloMEtPhi_[i] = -99.; 
00266                 ErsatzV1T1MEt_[i] = -99.; ErsatzV1T1Mt_[i] = -99.; ErsatzV1T1MEtPhi_[i] = -99.; 
00267                 ErsatzV1PfMEt_[i] = -99.; ErsatzV1PfMt_[i] = -99.; ErsatzV1PfMEtPhi_[i] = -99.; 
00268                 ErsatzV1TcMEt_[i] = -99.; ErsatzV1TcMt_[i] = -99.; ErsatzV1TcMEtPhi_[i] = -99.; 
00269                 
00270                 probe_sc_pt_[i] = -99.; probe_sc_eta_[i] = -99.; probe_sc_phi_[i] = -99.;
00271                 probe_sc_E_[i] = -99.; probe_sc_rawE_[i] = -99.;
00272                 probe_scV2_E_[i] = -99.;
00273                 probe_scV3_E_[i] = -99.;
00274                 probe_scV4_E_[i] = -99.;
00275 
00276                 if(Zevent_)
00277                 {
00278                         McElec_pt_[i] = -99.; McElec_eta_[i] = -99.; McElec_phi_[i] = -99.;
00279                         McElec_rescPt_[i] = -99.; McElec_rescEta_[i] = -99.; McElec_rescPhi_[i] = -99.;
00280                         McProbe_pt_[i] = -99.; McProbe_eta_[i] = -99.; McProbe_phi_[i] = -99.;
00281                         McProbe_rescPt_[i] = -99.; McProbe_rescEta_[i] = -99.; McProbe_rescPhi_[i] = -99.;
00282                         McElecProbe_dPhi_[i] = -99.; McElecProbe_dEta_[i] = -99.; McElecProbe_dR_[i] = -99.;
00283                 }
00284 
00285                 edm::LogDebug_("","",180)<<"Initialisation of array index "<< i <<" completed.";
00286         }       
00287         //Get Collections
00288         edm::Handle<reco::GenParticleCollection> pGenPart;
00289         try
00290         {
00291                 evt.getByLabel(MCTruthCollection_, pGenPart);
00292         }catch(cms::Exception &ex)
00293         {       
00294                 edm::LogError("analyze") <<"Can't get collection with label "<< MCTruthCollection_.label();
00295         } 
00296         edm::Handle<reco::GsfElectronCollection> pElectrons;
00297         try
00298         {
00299                 evt.getByLabel(ElectronCollection_, pElectrons);
00300         }catch(cms::Exception &ex)
00301         {
00302                 edm::LogError("analyze") <<"Can't get collection with label "<< ElectronCollection_.label();
00303         }
00304         edm::Handle<reco::SuperClusterCollection> pHybrid;
00305         try
00306         {
00307                 evt.getByLabel(HybridScCollection_, pHybrid);
00308         }catch(cms::Exception &ex)
00309         {
00310                 edm::LogError("analyze") <<"Can't get collection with label "<< HybridScCollection_.label();
00311         }
00312         edm::Handle<reco::SuperClusterCollection> pM5x5;
00313         try
00314         {
00315                 evt.getByLabel(M5x5ScCollection_, pM5x5);
00316         }catch(cms::Exception &ex)
00317         {
00318                 edm::LogError("analyze") <<"Can't get collection with label "<< M5x5ScCollection_.label();
00319         }
00320         edm::Handle<reco::CaloMETCollection> pCaloMEt;
00321         try
00322         {
00323                 evt.getByLabel(CaloMEtCollection_, pCaloMEt);
00324         }catch(cms::Exception &ex)
00325         {
00326                 edm::LogError("analyze")<<"Can't get collection with label "<< CaloMEtCollection_.label();
00327         }
00328         edm::Handle<reco::METCollection> pT1MEt;
00329 //      try
00330 //      {
00331 //              evt.getByLabel(T1MEtCollection_, pT1MEt);
00332 //      }catch(cms::Exception &ex)
00333 //      {
00334 //              edm::LogError("analyze")<<"Can't get collection with label "<< T1MEtCollection_.label();
00335 //      }
00336         edm::Handle<reco::PFMETCollection> pPfMEt;
00337         try
00338         {
00339                 evt.getByLabel(PfMEtCollection_, pPfMEt);
00340         }catch(cms::Exception &ex)
00341         {
00342                 edm::LogError("analyze")<<"Can't get collection with label "<< PfMEtCollection_.label();
00343         }
00344         edm::Handle<reco::METCollection> pTcMEt;
00345         try
00346         {
00347                 evt.getByLabel(TcMEtCollection_, pTcMEt);
00348         }catch(cms::Exception &ex)
00349         {
00350                 edm::LogError("analyze")<<"Can't get collection with label "<< TcMEtCollection_.label();
00351         }
00352         edm::Handle<reco::GenMETCollection> pGenMEt;
00353         try
00354         {
00355                 evt.getByLabel(GenMEtCollection_, pGenMEt);
00356         }catch(cms::Exception &ex)
00357         {
00358                 edm::LogError("analyze")<<"Can't get collection with label "<< GenMEtCollection_.label();
00359         }
00360         edm::Handle<edm::TriggerResults> pTriggerResults;
00361         try
00362         {
00363                 evt.getByLabel(TriggerResults_, pTriggerResults);
00364         }catch(cms::Exception &ex)
00365         {
00366                 edm::LogError("analyze")<<"Cant get collection with label "<< TriggerResults_.label();
00367         }
00368         edm::Handle<trigger::TriggerEvent> pHLT;
00369         try
00370         {
00371                 evt.getByLabel(TriggerEvent_, pHLT);
00372         }catch(cms::Exception &ex)
00373         {
00374                 edm::LogError("analyze")<<"Can't get collection with label "<< TriggerEvent_.label();
00375         }
00376 
00377         std::vector<math::XYZTLorentzVector>McElecs,McElecsFinalState;
00378         std::vector<math::XYZTLorentzVector> McElecsResc;
00379         if(Zevent_)
00380         {
00381                 edm::LogDebug_("","",289)<<"Analysing MC properties.";
00382                 const reco::GenParticleCollection *McCand = pGenPart.product();
00383                 math::XYZTLorentzVector Zboson, RescZboson, McElec1, McElec2; 
00384                 for(reco::GenParticleCollection::const_iterator McP = McCand->begin(); McP != McCand->end(); ++McP)
00385                 {
00386                         const reco::Candidate* mum = McP->mother();
00387                         if(abs(McP->pdgId())==11 && abs(mum->pdgId()) == 23)
00388                         {
00389                                 McElecs.push_back(McP->p4());
00390                                 if(abs(mum->pdgId() == 23)) Zboson = mum->p4();
00391                                 
00392                                 std::cout <<"Found electron, ID = "<< McP->pdgId() <<"\t status = "<< McP->status()<<std::endl;
00393                                 if(McP->status() != 1)
00394                                 {
00395                                         const reco::Candidate* McPD = McP->daughter(0);
00396                                         McPD = McPD->mother();
00397                                         while(McPD->status() != 1)
00398                                         {
00399                                                 int n = McPD->numberOfDaughters();
00400                                                 std::cout<< McPD->pdgId() << " : status = "<<McPD->status()
00401                                                         <<"\tNumber of Daughters = "<< n <<std::endl;
00402                                                 for(int j = 0; j < n; ++ j) 
00403                                                 {
00404                                                         const reco::Candidate *d = McPD->daughter( j );
00405                                                         std::cout <<"Daughter "<< j <<"\t id = "<< d->pdgId() << std::endl;
00406                                                         if(abs(d->pdgId()) == 11)
00407                                                         { 
00408                                                                 McPD = d;
00409                                                                 break;
00410                                                         }
00411                                                 }
00412                                         }
00413                                         std::cout<< McPD->pdgId() << " : status = "<<McPD->status()<<"\tAdding to vector!"<<std::endl;
00414                                         McElecsFinalState.push_back(McPD->p4());
00415                                 }else McElecsFinalState.push_back(McP->p4());
00416                         }
00417                 }
00418                 McZ_m_ = Zboson.M(); McZ_pt_ = Zboson.Pt(); McZ_phi_ = Zboson.Phi(); McZ_eta_ = Zboson.Eta(); McZ_y_ = Zboson.Rapidity();
00419                 McElec_nZmum_ =McElecs.size();
00420                 McElec_nFinal_ =McElecsFinalState.size();
00421                 edm::LogDebug_("","",309)<<"MC electrons with Z mother = "<< McElec_nZmum_ 
00422                                                 <<"\tFinal state MC electrons = "<< McElec_nFinal_;
00423                 
00424                 McElecsResc.resize(2);
00425 //              RescZboson.SetCoordinates(Zboson.Px(), Zboson.Py(), Zboson.Pz(), sqrt(Zboson.P2()+(mW_*mW_*Zboson.M2())/(mZ_*mZ_)));
00426                 RescZboson.SetCoordinates(Zboson.Px()*mW_/mZ_, Zboson.Py()*mW_/mZ_, Zboson.Pz()*mW_/mZ_, Zboson.E()*mW_/mZ_);
00427                 McZ_rescM_ = RescZboson.M(); McZ_rescPt_ = RescZboson.Pt(); McZ_rescEta_ = RescZboson.Eta(); McZ_rescPhi_ = RescZboson.Phi();
00428                 McZ_rescY_ = RescZboson.Rapidity();
00429                 ROOT::Math::Boost CoMBoost(Zboson.BoostToCM());
00430                 
00431                 math::XYZTLorentzVector RescMcElec0 = CoMBoost(McElecsFinalState[0]);
00432                 math::XYZTLorentzVector RescMcElec1 = CoMBoost(McElecsFinalState[1]);
00433                 RescMcElec0 *= mW_/mZ_;
00434                 RescMcElec1 *= mW_/mZ_;
00435 
00436                 double E_W = RescZboson.E();
00437                 ROOT::Math::Boost BackToLab(RescZboson.Px()/E_W, RescZboson.Py()/E_W, RescZboson.Pz()/E_W);
00438 
00439                 RescMcElec0 = BackToLab(RescMcElec0);
00440 //              RndmMcElec_Rescaled_pt_ = RescMcElec0.Pt();
00441 //              RndmMcElec_Rescaled_eta_ = RescMcElec0.Eta();
00442 //              RndmMcElec_Rescaled_phi_ = RescMcElec0.Phi();
00443         
00444                 RescMcElec1 = BackToLab(RescMcElec1);
00445 //              OthrMcElec_Rescaled_pt_ = RescMcElec1.Pt();
00446 //              OthrMcElec_Rescaled_eta_ = RescMcElec1.Eta();
00447 //              OthrMcElec_Rescaled_phi_ = RescMcElec1.Phi();
00448                 McElecsResc[0] = RescMcElec0;
00449                 McElecsResc[1] = RescMcElec1;
00450                 math::XYZTLorentzVector sum = RescMcElec1+RescMcElec0;
00451                 edm::LogDebug_("","", 307)<<"McElecsResc[0] + McElecsResc[1] = ("<<sum.Px()<<", "<<sum.Py()<<", "
00452                                                 <<sum.Pz()<<", "<<sum.E()<<")";
00453         }       
00454 
00455         const edm::TriggerResults* HltRes = pTriggerResults.product();
00456         const edm::TriggerNames & triggerNames = evt.triggerNames(*HltRes);
00457         if(HLTPathCheck_)
00458         {
00459                 for(unsigned int itrig = 0; itrig < HltRes->size(); ++itrig)
00460                 {
00461                         std::string nom = triggerNames.triggerName(itrig);
00462                         edm::LogInfo("")<< itrig <<" : Name = "<< nom <<"\t Accepted = "<< HltRes->accept(itrig);
00463                 }
00464         }
00465         if(HltRes->accept(34) ==0) edm::LogError("")<<"Event did not pass "<< triggerNames.triggerName(34)<<"!";
00466         if(HltRes->accept(34) !=0)
00467         {
00468         std::vector<reco::GsfElectronRef> UniqueElectrons;
00469         UniqueElectrons = uniqueElectronFinder(pElectrons);
00470         edm::LogDebug_("","ErsatzMEt",192)<<"Unique electron size = "<<UniqueElectrons.size();
00471         std::vector<reco::GsfElectronRef> SelectedElectrons;
00472         const unsigned int fId = pHLT->filterIndex(TriggerPath_);
00473         std::cout << "Filter Id = " << fId << std::endl;
00474         SelectedElectrons = electronSelector(UniqueElectrons, pHLT, fId, CutVector_);
00475         nTags_ = SelectedElectrons.size();
00476         edm::LogDebug_("","ErsatzMEt",197)<<"Selected electron size = "<<nTags_;
00477 
00478         iComb_ = 0;
00479         if(Zevent_)
00480         {
00481         //Match MC electrons to the selected electrons and store some of their properties in the tree.
00482         //The properties of the other MC electron (i.e. that not selected) are also stored.
00483                 for(std::vector<reco::GsfElectronRef>::const_iterator elec = SelectedElectrons.begin(); 
00484                                                                 elec != SelectedElectrons.end(); ++elec)
00485                 {
00486                         for(int m = 0; m < 2; ++m)
00487                         {
00488                                 double dRLimit = 99.;
00489                                 double dR = reco::deltaR(McElecs[m], *(*elec)); 
00490                                 if(dR < dRLimit)
00491                                 {
00492                                         dRLimit = dR;
00493                                         McElec_pt_[iComb_] = McElecs[m].pt();
00494                                         McElec_eta_[iComb_] = McElecs[m].eta();
00495                                         McElec_rescPt_[iComb_] = McElecsResc[m].pt();
00496                                         McElec_rescEta_[iComb_] = McElecsResc[m].eta();
00497                                 } 
00498                         }
00499                 }
00500         }
00501 
00502         std::map<reco::GsfElectronRef, reco::GsfElectronRef> TagProbePairs;
00503         TagProbePairs = probeFinder(SelectedElectrons, pElectrons);
00504         nProbes_ = TagProbePairs.size();        
00505         edm::LogDebug_("", "ErsatzMEt", 209)<<"Number of tag-probe pairs = "<< TagProbePairs.size();
00506 
00507         if(!TagProbePairs.empty())
00508         {
00509                 const reco::CaloMETCollection* caloMEtCollection = pCaloMEt.product();
00510                 const reco::MET calomet = *(caloMEtCollection->begin());
00511                 CaloMEt_ = calomet.pt();
00512                 CaloMEtphi_ = calomet.phi();
00513 
00514                 //const reco::METCollection* t1MEtCollection = pT1MEt.product();
00515                 //const reco::MET t1met = *(t1MEtCollection->begin());
00516                 //T1MEt_ = t1met.pt();
00517                 //T1MEtphi_ = t1met.phi();
00518                 
00519                 const reco::PFMETCollection* pfMEtCollection = pPfMEt.product();
00520                 const reco::MET pfmet = *(pfMEtCollection->begin());
00521                 PfMEt_ = pfmet.pt();
00522                 PfMEtphi_ = pfmet.phi();
00523                 
00524                 const reco::METCollection* tcMEtCollection = pTcMEt.product();
00525                 const reco::MET tcmet = *(tcMEtCollection->begin());
00526                 TcMEt_ = tcmet.pt();
00527                 TcMEtphi_ = tcmet.phi();
00528                 
00529                 reco::MET ersatzMEt;
00530 
00531                 for(std::map<reco::GsfElectronRef, reco::GsfElectronRef>::const_iterator it = TagProbePairs.begin();
00532                         it != TagProbePairs.end(); ++it)
00533                 {
00534                         edm::LogDebug_("","DelendumLoop", 293)<<"iComb_ = "<< iComb_;
00535                         tag_q_[iComb_] = it->first->charge();
00536                         edm::LogDebug_("","",360)<<"tag charge = "<< tag_q_[iComb_];
00537                         tag_pt_[iComb_] = it->first->pt();
00538                         tag_eta_[iComb_] = it->first->eta();
00539                         tag_phi_[iComb_] = it->first->phi();
00540                         edm::LogDebug_("","ErsatzMEt", 364)<<"tag pt = "<< tag_pt_[iComb_] 
00541                                         <<"\teta = "<< tag_eta_[iComb_]<<"\tphi = "<< tag_phi_[iComb_];
00542                         tag_trckIso_[iComb_] = it->first->isolationVariables03().tkSumPt;
00543                         tag_ecalIso_[iComb_] = it->first->isolationVariables04().ecalRecHitSumEt;
00544                         tag_hcalIso_[iComb_] = it->first->isolationVariables04().hcalDepth1TowerSumEt 
00545                                                 + it->first->isolationVariables04().hcalDepth2TowerSumEt;
00546                         edm::LogDebug_("","ErsatzMEt", 370)<<"tag trackiso = "<< tag_trckIso_[iComb_] 
00547                                         <<"\tecaliso = "<< tag_ecalIso_[iComb_]<<"\thcaliso = "<< tag_hcalIso_[iComb_];
00548                         tag_sIhIh_[iComb_] = it->first->scSigmaIEtaIEta();
00549                         tag_dPhiIn_[iComb_] = it->first->deltaPhiSuperClusterTrackAtVtx();
00550                         tag_dEtaIn_[iComb_] = it->first->deltaEtaSuperClusterTrackAtVtx();
00551                         edm::LogDebug_("","ErsatzMEt", 245)<<"tag sIhIh = "<< tag_sIhIh_[iComb_] 
00552                                         <<"\tdPhiIn = "<< tag_dPhiIn_[iComb_]<<"\tdEtaIn = "<< tag_dEtaIn_[iComb_];
00553                         tag_e5x5_[iComb_] = it->first->scE5x5();
00554                         tag_e2x5Max_[iComb_] = it->first->scE2x5Max();
00555                         tag_e2x5Max_[iComb_] = it->first->scE1x5();
00556                         edm::LogDebug_("","ErsatzMEt", 245)<<"tag e5x5 = "<< tag_e5x5_[iComb_] 
00557                                         <<"\te2x5Max = "<< tag_e2x5Max_[iComb_]<<"\te1x5Max = "<< tag_e1x5Max_[iComb_];
00558                         tag_hoe_[iComb_] = it->first->hadronicOverEm();
00559                         tag_eop_[iComb_] = it->first->eSuperClusterOverP();
00560                         tag_pin_[iComb_] = it->first->trackMomentumAtVtx().R();
00561                         tag_pout_[iComb_] = it->first->trackMomentumOut().R();
00562                         edm::LogDebug_("","ErsatzMEt", 245)<<"tag hoe = "<<tag_hoe_[iComb_]<<"\tpoe = "<<tag_eop_[iComb_] 
00563                                         <<"\tpin = "<< tag_pin_[iComb_]<<"\tpout = "<< tag_pout_[iComb_];
00564                         probe_q_[iComb_] = it->first->charge();
00565                         edm::LogDebug_("","",360)<<"probe charge = "<< probe_q_[iComb_];
00566                         probe_pt_[iComb_] = it->second->pt();
00567                         probe_eta_[iComb_] = it->second->eta();
00568                         probe_phi_[iComb_] = it->second->phi();
00569                         edm::LogDebug_("","ErsatzMEt", 245)<<"probe pt = "<< probe_pt_[iComb_] 
00570                                         <<"\teta = "<< probe_eta_[iComb_]<<"\tphi = "<< probe_phi_[iComb_];
00571                         probe_trckIso_[iComb_] = it->second->isolationVariables03().tkSumPt;
00572                         probe_ecalIso_[iComb_] = it->second->isolationVariables04().ecalRecHitSumEt;
00573                         probe_hcalIso_[iComb_] = it->second->isolationVariables04().hcalDepth1TowerSumEt 
00574                                                 + it->second->isolationVariables04().hcalDepth2TowerSumEt;
00575                         edm::LogDebug_("","ErsatzMEt", 245)<<"probe trackiso = "<< probe_trckIso_[iComb_] 
00576                                         <<"\tecaliso = "<< probe_ecalIso_[iComb_]<<"\thcaliso = "<< probe_phi_[iComb_];
00577                         probe_sIhIh_[iComb_] = it->second->scSigmaIEtaIEta();
00578                         probe_dPhiIn_[iComb_] = it->second->deltaPhiSuperClusterTrackAtVtx();
00579                         probe_dEtaIn_[iComb_] = it->second->deltaEtaSuperClusterTrackAtVtx();
00580                         edm::LogDebug_("","ErsatzMEt", 245)<<"probe sIhIh = "<< probe_sIhIh_[iComb_] 
00581                                         <<"\tdPhiIn = "<< probe_dPhiIn_[iComb_]<<"\tdEtaIn = "<< probe_dEtaIn_[iComb_];
00582                         probe_e5x5_[iComb_] = it->second->scE5x5();
00583                         probe_e2x5Max_[iComb_] = it->second->scE2x5Max();
00584                         probe_e2x5Max_[iComb_] = it->second->scE1x5();
00585                         edm::LogDebug_("","ErsatzMEt", 245)<<"probe e5x5 = "<< probe_e5x5_[iComb_] 
00586                                         <<"\te2x5Max = "<< probe_e2x5Max_[iComb_]<<"\te1x5Max = "<< probe_e1x5Max_[iComb_];
00587                         probe_hoe_[iComb_] = it->second->hadronicOverEm();
00588                         probe_eop_[iComb_] = it->second->eSuperClusterOverP();
00589                         probe_pin_[iComb_] = it->second->trackMomentumAtVtx().R();
00590                         probe_pout_[iComb_] = it->second->trackMomentumOut().R();
00591                         edm::LogDebug_("","ErsatzMEt", 245)<<"probe hoe = "<<probe_hoe_[iComb_]<<"\tpoe = "<<probe_eop_[iComb_] 
00592                                         <<"\tpin = "<< probe_pin_[iComb_]<<"\tpout = "<< probe_pout_[iComb_];
00593                         
00594                         double dRLimit = 0.2;
00595                         for(unsigned int mcEId = 0; mcEId < McElecs.size(); ++mcEId)
00596                         {
00597 //                              double dR = reco::deltaR((*(*mcEl)), probeVec); 
00598                                 double dR = reco::deltaR(McElecs[mcEId], it->second->p4()); 
00599                                 if(dR < dRLimit)
00600                                 {
00601                                         dRLimit = dR;
00602                                         McProbe_pt_[iComb_] = McElecs[mcEId].pt();
00603                                         McProbe_eta_[iComb_] = McElecs[mcEId].eta();
00604                                         McProbe_phi_[iComb_] = McElecs[mcEId].phi();
00605                                         McProbe_rescPt_[iComb_] = McElecsResc[mcEId].pt();
00606                                         McProbe_rescEta_[iComb_] = McElecsResc[mcEId].eta();
00607                                         McProbe_rescPhi_[iComb_] = McElecsResc[mcEId].phi();
00608                                         probe_d_MCE_SCE_[iComb_] = McElecs[mcEId].energy() - it->second->superCluster()->rawEnergy();
00609                                         McElecProbe_dPhi_[iComb_] = reco::deltaPhi(McElecs[mcEId].phi(), McElecs[(mcEId+1)%2].phi());
00610                                         McElecProbe_dEta_[iComb_] = fabs(McElecs[mcEId].eta() - McElecs[(mcEId+1)%2].eta());
00611                                         McElecProbe_dR_[iComb_] = reco::deltaR(McElecs[mcEId], McElecs[(mcEId+1)%2]);
00612                                 } 
00613                         }
00614 
00615                         // Uncorrected supercluster V1
00616                         reco::SuperCluster scV1 = *(it->second->superCluster());
00617                         math::XYZTLorentzVector probe_scV1_detVec = DetectorVector(scV1);
00618                         probe_sc_pt_[iComb_] = probe_scV1_detVec.pt(); 
00619                         probe_sc_eta_[iComb_] = scV1.eta();
00620                         probe_sc_phi_[iComb_] = scV1.phi();
00621                         probe_sc_nClus_[iComb_] = scV1.clustersSize();
00622                         probe_sc_E_[iComb_] = scV1.energy();
00623                         probe_sc_rawE_[iComb_] = scV1.rawEnergy();
00624 
00625                         ersatzMEt = ersatzFabrik(it->first, scV1, calomet, 1);
00626                         ErsatzV1CaloMEt_[iComb_] = ersatzMEt.pt();
00627                         ErsatzV1CaloMEtPhi_[iComb_] = ersatzMEt.phi();
00628                         //ersatzMEt = ersatzFabrik(it->first, it->second, t1met);
00629                         //ErsatzV1T1MEt_[iComb_] = ersatzMEt.pt();
00630                         //ErsatzV1T1MEtPhi_[iComb_] = ersatzMEt.phi();
00631                         ersatzMEt = ersatzFabrik(it->first, it->second, pfmet);
00632                         ErsatzV1PfMEt_[iComb_] = ersatzMEt.pt();
00633                         ErsatzV1PfMEtPhi_[iComb_] = ersatzMEt.phi();
00634                         ersatzMEt = ersatzFabrik(it->first, it->second, tcmet);
00635                         ErsatzV1TcMEt_[iComb_] = ersatzMEt.pt();
00636                         ErsatzV1TcMEtPhi_[iComb_] = ersatzMEt.phi();
00637 
00638                         // fEta corrected supercluster V2
00639                         reco::SuperCluster scV2;
00640                         if(fabs(probe_sc_eta_[iComb_]) < 1.479) 
00641                         {
00642                                 scV2 = fEtaScCorr(scV1);
00643                         }else{
00644                                 scV2 = scV1;
00645                         }
00646                         probe_scV2_E_[iComb_] = scV2.energy();
00647                         ersatzMEt = ersatzFabrik(it->first, scV2, calomet, 2);
00648                         ErsatzV2CaloMEt_[iComb_] = ersatzMEt.pt();
00649                         ErsatzV2CaloMEtPhi_[iComb_] = ersatzMEt.phi();
00650                 
00651                         // fBrem corrected supercluster V3
00652                         reco::SuperCluster scV3;        
00653                         if(fabs(probe_sc_eta_[iComb_]) < 1.479) 
00654                         {
00655                                 scV3 = fBremScCorr(scV1, hyb_fCorrPSet_);
00656                         }else{
00657                                 scV3 = fBremScCorr(scV1, m5x5_fCorrPSet_);
00658                         }       
00659                         probe_scV3_E_[iComb_] = scV3.energy();
00660                         ersatzMEt = ersatzFabrik(it->first, scV3, calomet, 3);
00661                         ErsatzV3CaloMEt_[iComb_] = ersatzMEt.pt();
00662                         ErsatzV3CaloMEtPhi_[iComb_] = ersatzMEt.phi();
00663                         
00664                         // Fully corrected supercluster V4
00665                         reco::SuperCluster scV4;
00666                         if(fabs(probe_sc_eta_[iComb_]) < 1.479) 
00667                         {
00668                                 scV4 = fBremScCorr(scV1, hyb_fCorrPSet_);
00669                         }else{
00670                                 scV4 = fBremScCorr(scV1, m5x5_fCorrPSet_);
00671                         }       
00672                         probe_scV4_E_[iComb_] = scV4.energy();
00673                         ersatzMEt = ersatzFabrik(it->first, scV4, calomet, 4);
00674                         ErsatzV4CaloMEt_[iComb_] = ersatzMEt.pt();
00675                         ErsatzV4CaloMEtPhi_[iComb_] = ersatzMEt.phi();
00676 
00677                         ++iComb_;
00678                 }
00679                 t_->Fill();
00680         }
00681         }
00682 }
00683 
00684 std::map<reco::GsfElectronRef, reco::GsfElectronRef> ErsatzMEt::probeFinder(const std::vector<reco::GsfElectronRef>& tags,
00685                                                         const edm::Handle<reco::GsfElectronCollection> pElectrons)
00686 {
00687         const reco::GsfElectronCollection *probeCands = pElectrons.product();
00688         std::map<reco::GsfElectronRef, reco::GsfElectronRef> TagProbes;
00689         for(std::vector<reco::GsfElectronRef>::const_iterator tagelec = tags.begin(); tagelec != tags.end(); ++tagelec)
00690         {
00691                 reco::GsfElectronRef tag = *tagelec;
00692                 std::pair<reco::GsfElectronRef, reco::GsfElectronRef> TagProbePair;
00693                 int nProbesPerTag = 0;
00694                 int index = 0;
00695                 for(reco::GsfElectronCollection::const_iterator probeelec = probeCands->begin(); probeelec != probeCands->end(); ++probeelec)
00696                 {
00697                         reco::GsfElectronRef probe(pElectrons, index);
00698                         double probeScEta = probe->superCluster()->eta();
00699                         if(probe->superCluster() != tag->superCluster() && fabs(probeScEta) < 2.5)
00700                         {
00701                                 if(fabs(probeScEta) < 1.4442 || fabs(probeScEta) > 1.560)
00702                                 {
00703                                         double invmass = ROOT::Math::VectorUtil::InvariantMass(tag->p4(), probe->p4());
00704                                         if(mTPmin_ <= invmass && invmass <= mTPmax_)
00705                                         {
00706                                                 TagProbePair = std::make_pair(tag, probe);
00707                                                 ++nProbesPerTag;
00708                                         }               
00709                                 }
00710                         }
00711                         ++index;
00712                 }
00713                 //nGsfElectrons_ = index;
00714                 if(nProbesPerTag == 1) TagProbes.insert(TagProbePair);
00715         }       
00716         return TagProbes;
00717 }
00718 
00719 reco::MET ErsatzMEt::ersatzFabrik(const reco::GsfElectronRef& elec, 
00720                                         const reco::SuperCluster& sc,
00721                                         const reco::MET& met,
00722                                         const int corr)
00723 {
00724         const math::XYZPoint ZVertex(elec->TrackPositionAtVtx().X(), elec->TrackPositionAtVtx().Y(),elec->TrackPositionAtVtx().Z());
00725 
00726         math::XYZTLorentzVector nu, boost_nu, ele, boost_ele;
00727         reco::SuperCluster elecSc = *(elec->superCluster());
00728         nu = PhysicsVectorRaw(met.vertex(), sc);
00729         boost_nu = PhysicsVectorRaw(ZVertex, sc);
00730         ele = PhysicsVectorRaw(met.vertex(), elecSc);
00731         boost_ele = ele;
00732 
00733         //Should use reco vertex for best Z->ee measurement. 
00734         edm::LogDebug_("ersatzFabrikV1", "", 569)<<"elec  = ("<< elec->p4().Px() << ", "<< elec->p4().Py()<< ", "<< elec->p4().Pz() << ", "<< elec->p4().E()<<")";
00735         math::XYZTLorentzVector Zboson = boost_nu + elec->p4();
00736         edm::LogDebug_("ersatzFabrikV1", "", 569)<<"Z pt = "<< Zboson.Pt() << "Z boson mass = " << Zboson.M();
00737         edm::LogDebug_("ersatzFabrikV1","", 570)<<"Z boson in lab frame = ("<<Zboson.Px()<<", "<<Zboson.Py()<<", "
00738                                         <<Zboson.Pz()<<", "<<Zboson.E()<<")";
00739         math::XYZTLorentzVector RescZboson(Zboson.Px(), Zboson.Py(), Zboson.Pz(), sqrt(Zboson.P2()+(mW_*mW_*Zboson.M2())/(mZ_*mZ_)));
00740         edm::LogDebug_("ersatzFabrikV1","", 573)<<"W boson in lab frame = ("<<RescZboson.Px()<<", "<<RescZboson.Py()<<", "
00741                                         <<RescZboson.Pz()<<", "<<RescZboson.E()<<")";
00742         ROOT::Math::Boost BoostToZRestFrame(Zboson.BoostToCM());
00743         edm::LogDebug_("ersatzFabrikV1","", 576)<<"Electron in lab frame = ("<< boost_ele.Px()<<", "<< boost_ele.Py()<<", "
00744                                         << boost_ele.Pz()<<", "<< boost_ele.E()<<")";
00745         edm::LogDebug_("ersatzFabrikV1","", 578)<<"Ersatz Neutrino in lab frame = ("<< boost_nu.Px()<<", "<< boost_nu.Py()<<", "
00746                                         << boost_nu.Pz()<<", "<< boost_nu.E()<<")";
00747         boost_ele = BoostToZRestFrame(boost_ele);
00748         boost_nu = BoostToZRestFrame(boost_nu);
00749         edm::LogDebug_("ersatzFabrikV1","", 582)<<"Electron in Z rest frame = ("<<boost_ele.Px()<<", "<<boost_ele.Py()<<", "
00750                                         <<boost_ele.Pz()<<", "<<boost_ele.E()<<")";
00751         edm::LogDebug_("ersatzFabrikV1","", 584)<<"Ersatz Neutrino in Z rest frame = ("<<boost_nu.Px()<<", "<<boost_nu.Py()<<", "
00752                                         <<boost_nu.Pz()<<", "<<boost_nu.E()<<")";
00753         boost_ele *= mW_/mZ_;
00754         boost_nu *= mW_/mZ_;
00755 
00756         double E_W = RescZboson.E();
00757         ROOT::Math::Boost BackToLab(RescZboson.Px()/E_W, RescZboson.Py()/E_W, RescZboson.Pz()/E_W);
00758         math::XYZTLorentzVector metVec(-99999., -99999., -99., -99999.);
00759         boost_ele = BackToLab(boost_ele);
00760 
00761         boost_nu = BackToLab(boost_nu);
00762         math::XYZTLorentzVector sum = boost_nu+boost_ele;
00763         edm::LogDebug_("ersatzFabrikV1","", 597)<<"Electron back in lab frame = ("<<boost_ele.Px()<<", "<<boost_ele.Py()<<", "
00764                                         <<boost_ele.Pz()<<", "<<boost_ele.E()<<")";
00765         edm::LogDebug_("ersatzFabrikV1","", 599)<<"Ersatz Neutrino back in lab frame = ("<<boost_nu.Px()<<", "<<boost_nu.Py()<<", "
00766                                        <<boost_nu.Pz()<<", "<<boost_nu.E()<<")";
00767         edm::LogDebug_("ersatzFabrikV1","", 601)<<"boost_ele + boost_nu = ("<<sum.Px()<<", "<<sum.Py()<<", "
00768                                         <<sum.Pz()<<", "<<sum.E()<<")";
00769 
00770         nu.SetXYZT(nu.X(), nu.Y(), 0., nu.T());
00771         ele.SetXYZT(ele.X(), ele.Y(), 0., ele.T());
00772         boost_ele.SetXYZT(boost_ele.X(), boost_ele.Y(), 0., boost_ele.T());
00773         metVec = met.p4() + nu + ele - boost_ele;
00774         
00775         reco::MET ersatzMEt(metVec, met.vertex());      
00776         if (corr == 1)  
00777         {
00778                 //Z_caloV1_m_[iComb_] = Zboson.M();
00779                 //Z_caloV1_pt_[iComb_] = Zboson.Pt();
00780                 //Z_caloV1_y_[iComb_] = Zboson.Y();
00781                 //Z_caloV1_eta_[iComb_] = Zboson.Eta();
00782                 //Z_caloV1_phi_[iComb_] = Zboson.Phi();
00783                 //Z_caloV1_rescM_[iComb_] = RescZboson.M();
00784                 //Z_caloV1_rescPt_[iComb_] = RescZboson.Pt();
00785                 //Z_caloV1_rescY_[iComb_] = RescZboson.Y();
00786                 //Z_caloV1_rescEta_[iComb_] = RescZboson.Eta();
00787                 //Z_caloV1_rescPhi_[iComb_] = RescZboson.Phi();
00788                 //Z_caloV1_probe_dPhi_[iComb_] = reco::deltaPhi(Zboson.Phi(), elec->phi());
00789                 //tag_caloV1_rescPt_[iComb_]  = boost_ele.Pt();
00790                 //tag_caloV1_rescEta_[iComb_]  = boost_ele.Eta();
00791                 //tag_caloV1_rescPhi_[iComb_]  = boost_ele.Phi();
00792                 //probe_caloV1_rescPt_[iComb_]  = boost_nu.Pt();
00793                 //probe_caloV1_rescEta_[iComb_]  = boost_nu.Eta();
00794                 //probe_caloV1_rescPhi_[iComb_]  = boost_nu.Phi();
00795                 ErsatzV1_Mesc_[iComb_]  = ROOT::Math::VectorUtil::InvariantMass(elec->p4(), boost_nu);
00796                 ErsatzV1_rescMesc_[iComb_] = ROOT::Math::VectorUtil::InvariantMass(ele, nu);
00797                 ErsatzV1CaloMt_[iComb_] = sqrt(2.*boost_ele.Pt()*ersatzMEt.pt()*
00798                                                 (1-cos(reco::deltaPhi(boost_ele.Phi(), ersatzMEt.phi()))));
00799         }
00800         if (corr == 2)  
00801         {
00802                 //Z_caloV2_m_[iComb_] = Zboson.M();
00803                 //Z_caloV2_pt_[iComb_] = Zboson.Pt();
00804                 //Z_caloV2_y_[iComb_] = Zboson.Y();
00805                 //Z_caloV2_eta_[iComb_] = Zboson.Eta();
00806                 //Z_caloV2_phi_[iComb_] = Zboson.Phi();
00807                 //Z_caloV2_rescM_[iComb_] = RescZboson.M();
00808                 //Z_caloV2_rescPt_[iComb_] = RescZboson.Pt();
00809                 //Z_caloV2_rescY_[iComb_] = RescZboson.Y();
00810                 //Z_caloV2_rescEta_[iComb_] = RescZboson.Eta();
00811                 //Z_caloV2_rescPhi_[iComb_] = RescZboson.Phi();
00812                 //Z_caloV2_probe_dPhi_[iComb_] = reco::deltaPhi(Zboson.Phi(), boost_elec->phi());
00813                 //tag_caloV2_rescPt_[iComb_]  = boost_ele.Pt();
00814                 //tag_caloV2_rescEta_[iComb_]  = boost_ele.Eta();
00815                 //tag_caloV2_rescPhi_[iComb_]  = boost_ele.Phi();
00816                 //probe_caloV2_rescPt_[iComb_]  = boost_nu.Pt();
00817                 //probe_caloV2_rescEta_[iComb_]  = boost_nu.Eta();
00818                 //probe_caloV2_rescPhi_[iComb_]  = boost_nu.Phi();
00819                 ErsatzV2_Mesc_[iComb_]  = ROOT::Math::VectorUtil::InvariantMass(elec->p4(), boost_nu);
00820                 ErsatzV2_rescMesc_[iComb_] = ROOT::Math::VectorUtil::InvariantMass(ele, nu);
00821                 ErsatzV2CaloMt_[iComb_] = sqrt(2.*boost_ele.Pt()*ersatzMEt.pt()*
00822                                                 (1-cos(reco::deltaPhi(boost_ele.Phi(), ersatzMEt.phi()))));
00823         }
00824         if (corr == 3)  
00825         {
00826                 //Z_caloV3_m_[iComb_] = Zboson.M();
00827                 //Z_caloV3_pt_[iComb_] = Zboson.Pt();
00828                 //Z_caloV3_y_[iComb_] = Zboson.Y();
00829                 //Z_caloV3_eta_[iComb_] = Zboson.Eta();
00830                 //Z_caloV3_phi_[iComb_] = Zboson.Phi();
00831                 //Z_caloV3_rescM_[iComb_] = RescZboson.M();
00832                 //Z_caloV3_rescPt_[iComb_] = RescZboson.Pt();
00833                 //Z_caloV3_rescY_[iComb_] = RescZboson.Y();
00834                 //Z_caloV3_rescEta_[iComb_] = RescZboson.Eta();
00835                 //Z_caloV3_rescPhi_[iComb_] = RescZboson.Phi();
00836                 //Z_caloV3_probe_dPhi_[iComb_] = reco::deltaPhi(Zboson.Phi(), boost_elec->phi());
00837                 //tag_caloV3_rescPt_[iComb_]  = boost_ele.Pt();
00838                 //tag_caloV3_rescEta_[iComb_]  = boost_ele.Eta();
00839                 //tag_caloV3_rescPhi_[iComb_]  = boost_ele.Phi();
00840                 //probe_caloV3_rescPt_[iComb_]  = boost_nu.Pt();
00841                 //probe_caloV3_rescEta_[iComb_]  = boost_nu.Eta();
00842                 //probe_caloV3_rescPhi_[iComb_]  = boost_nu.Phi();
00843                 ErsatzV3_Mesc_[iComb_]  = ROOT::Math::VectorUtil::InvariantMass(elec->p4(), boost_nu);
00844                 ErsatzV3_rescMesc_[iComb_] = ROOT::Math::VectorUtil::InvariantMass(ele, nu);
00845                 ErsatzV3CaloMt_[iComb_] = sqrt(2.*boost_ele.Pt()*ersatzMEt.pt()*
00846                                                 (1-cos(reco::deltaPhi(boost_ele.Phi(), ersatzMEt.phi()))));
00847         }
00848         if (corr == 4)  
00849         {
00850                 //Z_caloV4_m_[iComb_] = Zboson.M();
00851                 //Z_caloV4_pt_[iComb_] = Zboson.Pt();
00852                 //Z_caloV4_y_[iComb_] = Zboson.Y();
00853                 //Z_caloV4_eta_[iComb_] = Zboson.Eta();
00854                 //Z_caloV4_phi_[iComb_] = Zboson.Phi();
00855                 //Z_caloV4_rescM_[iComb_] = RescZboson.M();
00856                 //Z_caloV4_rescPt_[iComb_] = RescZboson.Pt();
00857                 //Z_caloV4_rescY_[iComb_] = RescZboson.Y();
00858                 //Z_caloV4_rescEta_[iComb_] = RescZboson.Eta();
00859                 //Z_caloV4_rescPhi_[iComb_] = RescZboson.Phi();
00860                 //Z_caloV4_probe_dPhi_[iComb_] = reco::deltaPhi(Zboson.Phi(), boost_elec->phi());
00861                 //tag_caloV4_rescPt_[iComb_]  = boost_ele.Pt();
00862                 //tag_caloV4_rescEta_[iComb_]  = boost_ele.Eta();
00863                 //tag_caloV4_rescPhi_[iComb_]  = boost_ele.Phi();
00864                 //probe_caloV4_rescPt_[iComb_]  = boost_nu.Pt();
00865                 //probe_caloV4_rescEta_[iComb_]  = boost_nu.Eta();
00866                 //probe_caloV4_rescPhi_[iComb_]  = boost_nu.Phi();
00867                 ErsatzV4_Mesc_[iComb_]  = ROOT::Math::VectorUtil::InvariantMass(elec->p4(), boost_nu);
00868                 ErsatzV4_rescMesc_[iComb_] = ROOT::Math::VectorUtil::InvariantMass(ele, nu);
00869                 ErsatzV4CaloMt_[iComb_] = sqrt(2.*boost_ele.Pt()*ersatzMEt.pt()*
00870                                                 (1-cos(reco::deltaPhi(boost_ele.Phi(), ersatzMEt.phi()))));
00871         }
00872         return ersatzMEt;
00873 }
00874 
00875 reco::MET ErsatzMEt::ersatzFabrik(const reco::GsfElectronRef& tag, 
00876                                         const reco::GsfElectronRef& probe,
00877                                         const reco::MET& met)
00878 {
00879         math::XYZTLorentzVector elec, nu, boost_elec, boost_nu;
00880         boost_elec = tag->p4();
00881         edm::LogDebug_("ersatzFabrikV1", "", 858)<<"boost_elec  = ("<< boost_elec.Px() << ", "<< boost_elec.Py()<< ", "<< boost_elec.Pz() << ", "<< boost_elec.E()<<")";
00882         boost_nu = probe->p4();
00883         edm::LogDebug_("ersatzFabrikV1", "", 860)<<"boost_nu  = ("<< boost_nu.Px() << ", "<< boost_nu.Py()<< ", "<< boost_nu.Pz() << ", "<< boost_nu.E()<<")";
00884         math::XYZTLorentzVector Zboson = boost_elec + boost_nu;
00885         edm::LogDebug_("ersatzFabrikV1", "", 862)<<"Zboson  = ("<< Zboson.Px() << ", "<< Zboson.Py()<< ", "<< Zboson.Pz() << ", "<< Zboson.E()<<")";
00886         math::XYZTLorentzVector RescZboson(Zboson.Px(), Zboson.Py(), Zboson.Pz(), sqrt(Zboson.P2()+(mW_*mW_*Zboson.M2())/(mZ_*mZ_)));
00887         edm::LogDebug_("ersatzFabrikV1", "", 864)<<"RescZboson  = ("<< RescZboson.Px() << ", "<< RescZboson.Py()<< ", "<< RescZboson.Pz() << ", "<< RescZboson.E()<<")";
00888         ROOT::Math::Boost BoostToZRestFrame(Zboson.BoostToCM());
00889         elec = BoostToZRestFrame(boost_elec);
00890         edm::LogDebug_("ersatzFabrikV1", "", 867)<<"boost_elec (in Z rest frame) = ("<< elec.Px() << ", "<< elec.Py()<< ", "<< elec.Pz() << ", "<< elec.E()<<")";
00891         nu = BoostToZRestFrame(boost_nu);
00892         edm::LogDebug_("ersatzFabrikV1", "", 869)<<"boost_nu (in Z rest frame) = ("<< nu.Px() << ", "<< nu.Py()<< ", "<< nu.Pz() << ", "<< nu.E()<<")";
00893         elec *= mW_/mZ_;
00894         edm::LogDebug_("ersatzFabrikV1", "", 871)<<"elec (in Z rest frame) = ("<< elec.Px() << ", "<< elec.Py()<< ", "<< elec.Pz() << ", "<< elec.E()<<")";
00895         nu *= mW_/mZ_;
00896         edm::LogDebug_("ersatzFabrikV1", "", 873)<<"nu (in Z rest frame) = ("<< nu.Px() << ", "<< nu.Py()<< ", "<< nu.Pz() << ", "<< nu.E()<<")";
00897         ROOT::Math::Boost BoostBackToLab(RescZboson.Px()/RescZboson.E(), RescZboson.Py()/RescZboson.E(), RescZboson.Pz()/RescZboson.E());
00898         math::XYZTLorentzVector metVec(-99999., -99999., -99., -99999.);
00899         elec = BoostBackToLab(elec);
00900         edm::LogDebug_("ersatzFabrikV1", "", 877)<<"elec = ("<< elec.Px() << ", "<< elec.Py()<< ", "<< elec.Pz() << ", "<< elec.E()<<")";
00901         nu = BoostBackToLab(nu);
00902         edm::LogDebug_("ersatzFabrikV1", "", 879)<<"nu = ("<< nu.Px() << ", "<< nu.Py()<< ", "<< nu.Pz() << ", "<< nu.E()<<")";
00903         Z_m_[iComb_] = Zboson.M();
00904         Z_pt_[iComb_] = Zboson.Pt();
00905         Z_y_[iComb_] = Zboson.Y();
00906         Z_eta_[iComb_] = Zboson.Eta();
00907         Z_phi_[iComb_] = Zboson.Phi();
00908         Z_rescM_[iComb_] = RescZboson.M();
00909         Z_rescPt_[iComb_] = RescZboson.Pt();
00910         Z_rescY_[iComb_] = RescZboson.Y();
00911         Z_rescEta_[iComb_] = RescZboson.Eta();
00912         Z_rescPhi_[iComb_] = RescZboson.Phi();
00913         Z_probe_dPhi_[iComb_] = reco::deltaPhi(Zboson.Phi(), boost_elec.phi());
00914         tag_rescPt_[iComb_]  = elec.Pt();
00915         tag_rescEta_[iComb_]  = elec.Eta();
00916         tag_rescPhi_[iComb_]  = elec.Phi();
00917         probe_rescPt_[iComb_]  = nu.Pt();
00918         probe_rescEta_[iComb_]  = nu.Eta();
00919         probe_rescPhi_[iComb_]  = nu.Phi();
00920         elec.SetXYZT(elec.X(), elec.Y(), 0., elec.T());
00921         nu.SetXYZT(nu.X(), nu.Y(), 0., nu.T());
00922         boost_elec.SetXYZT(boost_elec.X(), boost_elec.Y(), 0., boost_elec.T());
00923         metVec = met.p4() + nu + elec - boost_elec;
00924         reco::MET ersatzMEt(metVec, met.vertex());      
00925         return ersatzMEt;
00926 }
00927 
00928 bool ErsatzMEt::isInBarrel(double eta)
00929 {
00930         return (fabs(eta) < BarrelEtaMax_);
00931 }
00932 
00933 bool ErsatzMEt::isInEndCap(double eta)
00934 {
00935         return (fabs(eta) < EndCapEtaMax_ && fabs(eta) > EndCapEtaMin_);
00936 }
00937 
00938 bool ErsatzMEt::isInFiducial(double eta)
00939 {
00940         return isInBarrel(eta) || isInEndCap(eta);
00941 }
00942 
00943 // ------------ method called once each job just after ending the event loop  ------------
00944 void ErsatzMEt::endJob() {
00945 }
00946 //define this as a plug-in
00947 DEFINE_FWK_MODULE(ErsatzMEt);
00948