// -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/RivetAIDA.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Tools/ParticleIdUtils.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/VisibleFinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/UnstableFinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { class ATLAS_2012_I1204447 : public Analysis { public: /// Constructor ATLAS_2012_I1204447() : Analysis("ATLAS_2012_I1204447") { setNeedsCrossSection(true); } // the fiducial efficiencies model the effects of the ATLAS detector bool use_fiducial_lepton_efficiency; // list of signal regions and event counts per signal region std::vector Signal_Regions; std::map EventCountsPerSR; /// Book histograms and initialise projections before the run void init() { // to calculate the acceptance without having the fiducial lepton efficiencies included, this part can be turned off use_fiducial_lepton_efficiency = true; // random numbers for simulation of ATLAS detector reconstruction efficiency srand (160385); // Read in all signal regions Signal_Regions = Get_Signal_Regions(); // set number of events per signal region to 0 for (unsigned int i=0; i < Signal_Regions.size(); i++) EventCountsPerSR[Signal_Regions.at(i)] = 0.0; /// final state including all charged and neutral particles const FinalState fs(-5.0, 5.0, 1000*MeV); addProjection(fs, "FS"); /// final state including all charged particles addProjection(ChargedFinalState(-2.5, 2.5, 1000*MeV), "CFS"); /// final state including all visible particles (to calculate MET, Jets etc.) addProjection(VisibleFinalState(-5.0,5.0),"vfs"); /// final state including all AntiKt 04 Jets VetoedFinalState vfs; vfs.addVetoPairId(MUON); addProjection(FastJets(vfs, FastJets::ANTIKT, 0.4), "AntiKtJets04"); /// final state including all unstable particles (including taus) addProjection(UnstableFinalState(-5.0, 5.0, 5.0*GeV),"UFS"); /// final state including all electrons IdentifiedFinalState elecs(-2.47,2.47, 10.0*GeV); elecs.acceptIdPair(ELECTRON); addProjection(elecs, "elecs"); /// final state including all muons IdentifiedFinalState muons(-2.5,2.5, 10.0*GeV); muons.acceptIdPair(MUON); addProjection(muons, "muons"); /// Book histograms: _h_HTlep_all = bookHistogram1D("HTlep_all" , 30, 0, 1500); _h_HTjets_all = bookHistogram1D("HTjets_all", 30, 0, 1500); _h_MET_all = bookHistogram1D("MET_all" , 20, 0, 1000); _h_Meff_all = bookHistogram1D("Meff_all" , 30, 0, 3000); _h_e_n = bookHistogram1D("e_n" , 10, -0.5, 9.5); _h_mu_n = bookHistogram1D("mu_n" , 10, -0.5, 9.5); _h_tau_n = bookHistogram1D("tau_n", 10, -0.5, 9.5); _h_pt_1_3l = bookHistogram1D("pt_1_3l", 100, 0, 2000); _h_pt_2_3l = bookHistogram1D("pt_2_3l", 100, 0, 2000); _h_pt_3_3l = bookHistogram1D("pt_3_3l", 100, 0, 2000); _h_pt_1_2ltau = bookHistogram1D("pt_1_2ltau", 100, 0, 2000); _h_pt_2_2ltau = bookHistogram1D("pt_2_2ltau", 100, 0, 2000); _h_pt_3_2ltau = bookHistogram1D("pt_3_2ltau", 100, 0, 2000); _h_excluded = bookHistogram1D("excluded", 2, -0.5, 1.5); } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = event.weight(); //---------- /// muons //---------- ParticleVector muon_candidate; ParticleVector charged_tracks = applyProjection(event, "CFS").particles(); ParticleVector visible_particles = applyProjection(event, "vfs").particles(); // loop over all muons foreach ( const Particle & mu, applyProjection(event, "muons").particlesByPt() ) { // calculate pTCone30 variable (pT of all tracks within dR<0.3 - pT of muon itself) double pTinCone = -mu.momentum().pT(); foreach ( const Particle & track, charged_tracks ) { if ( deltaR( mu.momentum(), track.momentum()) <= 0.3 ) pTinCone += track.momentum().pT(); } // calculate eTCone30 variable (pT of all visible particles within dR<0.3) double eTinCone = 0.; foreach ( const Particle & visible_particle, visible_particles ) { if ( deltaR( mu.momentum(), visible_particle.momentum()) < 0.3 && deltaR( mu.momentum(), visible_particle.momentum()) > 0.1 && abs( visible_particle.pdgId()) != 13 ) eTinCone += visible_particle.momentum().pT(); } /// apply reconstruction efficiency int muon_id = 13; if ( mu.hasAncestor(15) || mu.hasAncestor(-15)) muon_id = 14; // getting fiducial lepton efficiency double eff = 1.0; if (use_fiducial_lepton_efficiency) eff = apply_reco_eff( muon_id,mu.momentum().pT()/GeV, mu.momentum().eta()); // throwing dice if detector reconstructs particle bool keep_muon = rand()/static_cast(RAND_MAX) <= eff; // keeping muon if pTCone30/pT < 0.15 and eTCone30/pT < 0.2 and reconstructed by detector if ( pTinCone/mu.momentum().pT() <= 0.15 && eTinCone/mu.momentum().pT() < 0.2 && keep_muon ) muon_candidate.push_back(mu); } //----------- /// electrons //----------- ParticleVector electron_candidates; // loop over all electrons foreach ( const Particle & e, applyProjection(event, "elecs").particlesByPt() ) { // neglecting events in certain eta regions (crack regions) if( fabs(e.momentum().eta()) > 1.37 && fabs(e.momentum().eta()) < 1.52) continue; // calculate pTCone30 variable (pT of all tracks within dR<0.3 - pT of electron itself) double pTinCone = -e.momentum().perp(); foreach ( const Particle & track, charged_tracks ) { if ( deltaR( e.momentum(), track.momentum()) <= 0.3 ) pTinCone += track.momentum().pT(); } // calculate eTCone30 variable (pT of all visible particles (except muons) within dR<0.3) double eTinCone = 0.; foreach ( const Particle & visible_particle, visible_particles ) { if ( deltaR( e.momentum(), visible_particle.momentum()) < 0.3 && deltaR( e.momentum(), visible_particle.momentum()) > 0.1 && abs( visible_particle.pdgId() ) != 13 ) eTinCone += visible_particle.momentum().pT(); } /// apply reconstruction efficiency int elec_id = 11; if (e.hasAncestor(15) || e.hasAncestor(-15)) elec_id = 12; // getting fiducial lepton efficiency double eff = 1.0; if (use_fiducial_lepton_efficiency) eff = apply_reco_eff( elec_id, e.momentum().pT()/GeV, e.momentum().eta() ); // throwing dice if detector reconstructs particle bool keep_elec = rand()/static_cast(RAND_MAX) <= eff; // keeping electron if pTCone30/pT < 0.13 and eTCone30/pT < 0.2 and reconstructed by detector if ( pTinCone/e.momentum().perp() <= 0.13 && eTinCone/e.momentum().perp() < 0.2 && keep_elec) { electron_candidates.push_back(e); } } //------- /// taus //------- ParticleVector tau_candidates; // loop over all taus foreach ( const Particle & tau, applyProjection(event, "UFS").particlesByPt() ) { // only pick taus out of all unstable particles const int pid = tau.pdgId(); if ( abs(pid) != 15) continue; // get vertex of tau decaying into daughter particles const GenVertex* dv = tau.genParticle().end_vertex(); if (dv == 0) continue; /// calculate visible tau pT // get pT of tau neutrino in tau decay FourVector daughter_tau_neutrino_pT = get_tau_neutrino_pT(tau.genParticle()); // keep only taus in certain eta region and above 15 GeV of visible tau pT if ( (tau.genParticle().momentum() - daughter_tau_neutrino_pT).perp()/GeV <= 15.0 || fabs((tau.genParticle().momentum() - daughter_tau_neutrino_pT).eta()) > 2.5 ) continue; /// get prong number (number of tracks) in tau decay and check if tau decays leptonically unsigned int nprong = 0; bool lep_decaying_tau = false; get_prong_number(tau.genParticle(),nprong,lep_decaying_tau,tau.genParticle().barcode()); /// apply reconstruction efficiency int tau_id = 15; if (nprong == 1) tau_id = 15; if (nprong == 3) tau_id = 16; // getting fiducial lepton efficiency double eff = 1.0; if (use_fiducial_lepton_efficiency) eff = apply_reco_eff( tau_id, tau.momentum().pT()/GeV, tau.momentum().eta()); // throwing dice if detector reconstructs particle bool keep_tau = rand()/static_cast(RAND_MAX) <= eff; // keeping tau if nprong = 1 and tau decaying hadronically and reconstructed by detector if ( !lep_decaying_tau && nprong == 1 && keep_tau) tau_candidates.push_back(tau); } //------- /// jets //------- Jets jet_candidates; // getting all AntiKt04 Jets with pT > 25 GeV and eta < 4.9 foreach (const Jet& jet, applyProjection(event, "AntiKtJets04").jetsByPt(25.0*GeV) ) { if ( fabs( jet.momentum().eta() ) < 4.9 ) { jet_candidates.push_back(jet); } } //--------- /// ETmiss //--------- ParticleVector vfs_particles = applyProjection(event, "vfs").particles(); FourMomentum pTmiss; // loop over all visible particles, sum up the pT foreach (const Particle & p, vfs_particles ) { pTmiss -= p.momentum(); } // the difference to 0 is the ETmiss/MET value double eTmiss = pTmiss.pT()/GeV; //------------------ /// overlap removal //------------------ /// electron - electron ParticleVector electron_candidates_2; for (unsigned int ie=0; ie < electron_candidates.size(); ++ie) { const Particle & e = electron_candidates[ie]; bool away = true; // if electron pair within dR < 0.1: remove electron with lower pT for (unsigned int ie2=0; ie2 < electron_candidates_2.size(); ++ie2) { if ( deltaR( e.momentum(), electron_candidates_2[ie2].momentum()) < 0.1 ) { away = false; break; } } // if isolated keep it if ( away ) electron_candidates_2.push_back( e ); } /// jet - electron Jets recon_jets; foreach ( const Jet& jet, jet_candidates ) { bool away = true; // if jet within dR < 0.2 of electron: remove jet foreach ( const Particle & e, electron_candidates_2 ) { if ( deltaR( e.momentum(), jet.momentum()) < 0.2 ) { away = false; break; } } /// jet - tau if ( away ) { // if jet within dR < 0.2 of tau: remove jet foreach ( const Particle & tau, tau_candidates ) { if ( deltaR( tau.momentum(), jet.momentum()) < 0.2 ) { away = false; break; } } } // if isolated keep it if ( away ) recon_jets.push_back( jet ); } /// electron - jet ParticleVector recon_leptons; ParticleVector recon_e; for(unsigned int ie=0; ie < electron_candidates_2.size(); ++ie) { const Particle & e = electron_candidates_2[ie]; // if electron within 0.2 < dR < 0.4 from any jets: remove electron bool away = true; foreach ( const Jet& jet, recon_jets ) { if ( deltaR( e.momentum(), jet.momentum()) < 0.4 ) { away = false; break; } } /// electron - muon // if electron within dR < 0.1 of a muon: remove electron if ( away ) { foreach ( const Particle & mu, muon_candidate ) { if ( deltaR( mu.momentum(), e.momentum()) < 0.1 ) { away = false; break; } } } // if isolated keep it if ( away ) { recon_e.push_back( e ); recon_leptons.push_back( e ); } } /// tau - electron ParticleVector recon_tau; foreach ( const Particle & tau, tau_candidates ) { bool away = true; // if tau within dR < 0.2 of an electron: remove tau foreach ( const Particle & e, recon_e ) { if ( deltaR( tau.momentum(), e.momentum()) < 0.2 ) { away = false; break; } } /// tau - muon // if tau within dR < 0.2 of a muon: remove tau if ( away ) { foreach ( const Particle & mu, muon_candidate ) { if ( deltaR( tau.momentum(), mu.momentum()) < 0.2 ) { away = false; break; } } } // if isolated keep it if ( away ) recon_tau.push_back( tau ); } /// muon - jet ParticleVector recon_mu; ParticleVector trigger_mu; // if muon within dR < 0.4 of a jet: remove muon foreach ( const Particle & mu, muon_candidate ) { bool away = true; foreach ( const Jet& jet, recon_jets ) { if ( deltaR( mu.momentum(), jet.momentum()) < 0.4 ) { away = false; break; } } if ( away ) { recon_mu.push_back( mu ); recon_leptons.push_back( mu ); } if ( away && fabs(mu.momentum().eta())< 2.4) trigger_mu.push_back( mu ); } /// end overlap removal //--------------- /// jet cleaning //--------------- if( rand()/static_cast(RAND_MAX) <= 0.42) { foreach ( const Jet & jet, recon_jets ) { double eta = jet.momentum().rapidity(); double phi = jet.momentum().azimuthalAngle(MINUSPI_PLUSPI); if( jet.momentum().perp()/GeV > 25 && eta > -0.1 && eta < 1.5 && phi > -0.9 && phi < -0.5) vetoEvent; } } /// at least 3 charged tracks in event if (charged_tracks.size() < 3) vetoEvent; /// at least one e/mu passing trigger if( !( !recon_e .empty() && recon_e[0] .momentum().perp()/GeV > 25.) && !( !trigger_mu.empty() && trigger_mu[0].momentum().perp()/GeV > 25.) ) { MSG_DEBUG("Hardest lepton fails trigger"); vetoEvent; } /// only events with at least 2 electrons and muons and at least 3 leptons in total if ( recon_mu.size() + recon_e.size() + recon_tau.size() < 3 || recon_leptons.size() < 2) vetoEvent; std::sort(recon_leptons.begin(), recon_leptons.end(), cmpParticleByPt); std::sort(recon_tau.begin(), recon_tau.end(), cmpParticleByPt); /// calculate HTlep, fill lepton pT histograms & store chosen combination of 3 leptons double HTlep = 0.; ParticleVector chosen_leptons; if ( recon_leptons.size() > 2 ) { _h_pt_1_3l->fill(recon_leptons[0].momentum().perp()/GeV, weight); _h_pt_2_3l->fill(recon_leptons[1].momentum().perp()/GeV, weight); _h_pt_3_3l->fill(recon_leptons[2].momentum().perp()/GeV, weight); HTlep = recon_leptons[0].momentum().perp()/GeV + recon_leptons[1].momentum().perp()/GeV + recon_leptons[2].momentum().perp()/GeV; chosen_leptons.push_back( recon_leptons[0] ); chosen_leptons.push_back( recon_leptons[1] ); chosen_leptons.push_back( recon_leptons[2] ); } else { _h_pt_1_2ltau->fill(recon_leptons[0].momentum().perp()/GeV, weight); _h_pt_2_2ltau->fill(recon_leptons[1].momentum().perp()/GeV, weight); _h_pt_3_2ltau->fill(recon_tau[0].momentum().perp()/GeV, weight); HTlep = recon_leptons[0].momentum().perp()/GeV + recon_leptons[1].momentum().perp()/GeV + recon_tau[0].momentum().perp()/GeV ; chosen_leptons.push_back( recon_leptons[0] ); chosen_leptons.push_back( recon_leptons[1] ); chosen_leptons.push_back( recon_tau[0] ); } /// number of prompt e/mu and had taus _h_e_n ->fill(recon_e.size() , weight); _h_mu_n ->fill(recon_mu.size() , weight); _h_tau_n->fill(recon_tau.size(), weight); /// calculate HTjets double HTjets = 0.; foreach ( const Jet & jet, recon_jets ) HTjets += jet.momentum().perp()/GeV; /// calculate meff double meff = eTmiss + HTjets; ParticleVector all_leptons; foreach ( const Particle & e , recon_e ) { meff += e.momentum().perp()/GeV; all_leptons.push_back( e ); } foreach ( const Particle & mu, recon_mu ) { meff += mu.momentum().perp()/GeV; all_leptons.push_back( mu ); } foreach ( const Particle & tau, recon_tau ) { meff += tau.momentum().perp()/GeV; all_leptons.push_back( tau ); } /// fill histogram of kinematic variables _h_HTlep_all ->fill(HTlep , weight); _h_HTjets_all->fill(HTjets, weight); _h_MET_all ->fill(eTmiss, weight); _h_Meff_all ->fill(meff , weight); /// determine signal region (3l/2ltau, onZ/offZ) std::string basic_signal_region; if ( recon_mu.size() + recon_e.size() > 2 ) basic_signal_region += "3l_"; else if ( (recon_mu.size() + recon_e.size() == 2) && (recon_tau.size() > 0)) basic_signal_region += "2ltau_"; // is there an OSSF pair or a three lepton combination with an invariant mass close to the Z mass int onZ = isonZ(chosen_leptons); if (onZ == 1) basic_signal_region += "onZ"; else if (onZ == 0) basic_signal_region += "offZ"; // Check in which signal regions this event falls and adjust event counters Fill_EventCountsPerSR( basic_signal_region, onZ, HTlep, eTmiss, HTjets, meff, weight); } /// Normalise histograms etc., after the run void finalize() { // normalize to an integrated luminosity of 1 fb-1 double norm = crossSection()/femtobarn/sumOfWeights(); std::string best_signal_region = ""; double ratio_best_SR = 0.; // loop over all signal regions and find signal region with best sensitivity (ratio signal events/visible cross-section) for (unsigned int i=0; i < Signal_Regions.size(); i++) { double signal_events = EventCountsPerSR[Signal_Regions.at(i)] * norm; // use expected upper limits to find best signal region: double UL95 = Get_Upper_Limit( Signal_Regions.at(i), false); double ratio = signal_events / UL95; if (ratio > ratio_best_SR) { best_signal_region = Signal_Regions.at(i); ratio_best_SR = ratio; } } double signal_events_best_SR = EventCountsPerSR[best_signal_region] * norm; double exp_UL_best_SR = Get_Upper_Limit(best_signal_region, false); double obs_UL_best_SR = Get_Upper_Limit(best_signal_region, true); // print out result cout << "----------------------------------------------------------------------------------------" << endl; cout << "Best signal region: " << best_signal_region << endl; cout << "Normalized number of signal events in this best signal region (per fb-1): " << signal_events_best_SR << endl; cout << "Efficiency*Acceptance: " << EventCountsPerSR[best_signal_region]/sumOfWeights() << endl; cout << "Cross-section [fb]: " << crossSection()/femtobarn << endl; cout << "Expected visible cross-section (per fb-1): " << exp_UL_best_SR << endl; cout << "Ratio (signal events / expected visible cross-section): " << ratio_best_SR << endl; cout << "Observed visible cross-section (per fb-1): " << obs_UL_best_SR << endl; cout << "Ratio (signal events / observed visible cross-section): " << signal_events_best_SR/obs_UL_best_SR << endl; cout << "----------------------------------------------------------------------------------------" << endl; cout << "Using the EXPECTED limits (visible cross-section) of the analysis: " << endl; if (signal_events_best_SR > exp_UL_best_SR) { cout << "Since the number of signal events > the visible cross-section, this model/grid point is EXCLUDED with 95% CL." << endl; _h_excluded->fill(1); } else { cout << "Since the number of signal events < the visible cross-section, this model/grid point is NOT EXCLUDED." << endl; _h_excluded->fill(0); } cout << "----------------------------------------------------------------------------------------" << endl; cout << "Using the OBSERVED limits (visible cross-section) of the analysis: " << endl; if (signal_events_best_SR > obs_UL_best_SR) { cout << "Since the number of signal events > the visible cross-section, this model/grid point is EXCLUDED with 95% CL." << endl; _h_excluded->fill(1); } else { cout << "Since the number of signal events < the visible cross-section, this model/grid point is NOT EXCLUDED." << endl; _h_excluded->fill(0); } cout << "----------------------------------------------------------------------------------------" << endl; /// Normalize to cross section if (norm != 0) { scale(_h_HTlep_all, norm); scale(_h_HTjets_all, norm); scale(_h_MET_all, norm); scale(_h_Meff_all, norm); scale(_h_pt_1_3l, norm); scale(_h_pt_2_3l, norm); scale(_h_pt_3_3l, norm); scale(_h_pt_1_2ltau, norm); scale(_h_pt_2_2ltau, norm); scale(_h_pt_3_2ltau, norm); scale(_h_e_n, norm); scale(_h_mu_n, norm); scale(_h_tau_n, norm); scale(_h_excluded, signal_events_best_SR); } } private: /// Histograms AIDA::IHistogram1D *_h_HTlep_all; AIDA::IHistogram1D *_h_HTjets_all; AIDA::IHistogram1D *_h_MET_all; AIDA::IHistogram1D *_h_Meff_all; AIDA::IHistogram1D *_h_pt_1_3l; AIDA::IHistogram1D *_h_pt_2_3l; AIDA::IHistogram1D *_h_pt_3_3l; AIDA::IHistogram1D *_h_pt_1_2ltau; AIDA::IHistogram1D *_h_pt_2_2ltau; AIDA::IHistogram1D *_h_pt_3_2ltau; AIDA::IHistogram1D *_h_e_n; AIDA::IHistogram1D *_h_mu_n; AIDA::IHistogram1D *_h_tau_n; AIDA::IHistogram1D *_h_excluded; /// functions /// function giving a list of all signal regions std::vector Get_Signal_Regions() { // list of basic signal regions std::vector basic_signal_regions; basic_signal_regions.push_back("3l_offZ"); basic_signal_regions.push_back("3l_onZ"); basic_signal_regions.push_back("2ltau_offZ"); basic_signal_regions.push_back("2ltau_onZ"); // list of kinematic variables std::vector kinematic_variables; kinematic_variables.push_back("HTlep"); kinematic_variables.push_back("METStrong"); kinematic_variables.push_back("METWeak"); kinematic_variables.push_back("Meff"); kinematic_variables.push_back("MeffStrong"); std::vector signal_regions; // loop over all kinematic variables and basic signal regions for (unsigned int i0=0; i0 < kinematic_variables.size(); i0++) { for (unsigned int i1=0; i1 < basic_signal_regions.size(); i1++) { // is signal region onZ? int onZ = 0; std::size_t found = basic_signal_regions.at(i1).find("onZ"); if (found != std::string::npos) onZ = 1; // get cut values for this kinematic variable std::vector cut_values = Get_Cuts_Per_Signal_Region(kinematic_variables.at(i0), onZ); // loop over all cut values for (unsigned int i2=0; i2 < cut_values.size(); i2++) { // convert int to string std::ostringstream convert; convert << cut_values.at(i2); std::string cut_value = convert.str(); // push signal region into vector signal_regions.push_back( (kinematic_variables.at(i0) + "_" + basic_signal_regions.at(i1) + "_cut_" + cut_value).c_str() ); } } } return signal_regions; } /// function giving all cut vales per kinematic variable (taking onZ for MET into account) std::vector Get_Cuts_Per_Signal_Region(std::string signal_region, int onZ = 0) { std::vector Cut_Values; // cut values for HTlep if (signal_region.compare("HTlep") == 0) { Cut_Values.push_back(0); Cut_Values.push_back(100); Cut_Values.push_back(150); Cut_Values.push_back(200); Cut_Values.push_back(300); } // cut values for METStrong (HTjets > 100 GeV) and METWeak (HTjets < 100 GeV) else if (signal_region.compare("METStrong") == 0 || signal_region.compare("METWeak") == 0) { if (onZ == 0) Cut_Values.push_back(0); else if (onZ == 1) Cut_Values.push_back(20); Cut_Values.push_back(50); Cut_Values.push_back(75); } // cut values for Meff and MeffStrong (MET > 75 GeV) if (signal_region.compare("Meff") == 0 || signal_region.compare("MeffStrong") == 0) { Cut_Values.push_back(0); Cut_Values.push_back(150); Cut_Values.push_back(300); Cut_Values.push_back(500); } return Cut_Values; } /// function fills map EventCountsPerSR by looping over all signal regions /// and looking if the event falls into this signal region void Fill_EventCountsPerSR(std::string basic_signal_region, int onZ, double HTlep, double eTmiss, double HTjets, double meff, double weight) { // get cut values for HTlep, loop over them and add event if cut is passed std::vector cut_values = Get_Cuts_Per_Signal_Region("HTlep", onZ); for (unsigned int i=0; i < cut_values.size(); i++) { if (HTlep > cut_values.at(i)) EventCountsPerSR[("HTlep_" + basic_signal_region + "_cut_" + int2str(cut_values.at(i))).c_str()] += weight; } // get cut values for METStrong, loop over them and add event if cut is passed cut_values = Get_Cuts_Per_Signal_Region("METStrong", onZ); for (unsigned int i=0; i < cut_values.size(); i++) { if (eTmiss > cut_values.at(i) && HTjets > 100.) EventCountsPerSR[("METStrong_" + basic_signal_region + "_cut_" + int2str(cut_values.at(i))).c_str()] += weight; } // get cut values for METWeak, loop over them and add event if cut is passed cut_values = Get_Cuts_Per_Signal_Region("METWeak", onZ); for (unsigned int i=0; i < cut_values.size(); i++) { if (eTmiss > cut_values.at(i) && HTjets <= 100.) EventCountsPerSR[("METWeak_" + basic_signal_region + "_cut_" + int2str(cut_values.at(i))).c_str()] += weight; } // get cut values for Meff, loop over them and add event if cut is passed cut_values = Get_Cuts_Per_Signal_Region("Meff", onZ); for (unsigned int i=0; i < cut_values.size(); i++) { if (meff > cut_values.at(i)) EventCountsPerSR[("Meff_" + basic_signal_region + "_cut_" + int2str(cut_values.at(i))).c_str()] += weight; } // get cut values for MeffStrong, loop over them and add event if cut is passed cut_values = Get_Cuts_Per_Signal_Region("MeffStrong",onZ); for (unsigned int i=0; i < cut_values.size(); i++) { if (meff > cut_values.at(i) && eTmiss > 75.) EventCountsPerSR[("MeffStrong_" + basic_signal_region + "_cut_" + int2str(cut_values.at(i))).c_str()] += weight; } } /// function converting int to string std::string int2str (int Integer) { std::ostringstream convert; convert << Integer; std::string String = convert.str(); return String; } /// function returning 4-vector of daughter-particle if it is a tau neutrino FourVector get_tau_neutrino_pT(const GenParticle & p) { FourVector tau_neutrino_pT; const GenVertex* dv = p.end_vertex(); if (dv == 0) return p.momentum(); // loop over all daughter particles for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin(); pp != dv->particles_out_const_end(); ++pp) { int id = (*pp)->pdg_id(); if ( abs(id) == 16 ) tau_neutrino_pT = (*pp)->momentum(); } return tau_neutrino_pT; } /// function calculating the prong number of taus void get_prong_number(const GenParticle & p, unsigned int& nprong, bool& lep_decaying_tau, const int tau_barcode) { const GenVertex* dv = p.end_vertex(); if (dv != 0) { // loop over all daughter particles for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin(); pp != dv->particles_out_const_end(); ++pp) { int id = (*pp)->pdg_id(); // if they have status 1 and are charged they will produce a track and the prong number is +1 if ((*pp)->status() == 1 ) { if (Rivet::PID::charge(id) != 0 ) ++nprong; // check if tau decays leptonically if (( abs(id) == 11 || abs(id) == 13 || abs(id) == 15 ) && abs(p.pdg_id()) == 15) lep_decaying_tau = true; } // if the status of the daughter particle is 2 it is unstable and the further decays are checked else if ((*pp)->status() == 2 ) { get_prong_number((**pp),nprong,lep_decaying_tau, tau_barcode); } } } } /// function giving fiducial lepton efficiency double apply_reco_eff(int flavor, float pt, float eta) { double eff = 0.; double err = 0.; if (flavor == 11) { // weight prompt electron -- now including data/MC ID SF in eff. float rho = 0.820; float p0 = 7.34; float p1 = 0.8977; float ep0= 0.5 ; float ep1= 0.0087; eff = p1 - p0/pt; double err0 = ep0/pt; // d(eff)/dp0 double err1 = ep1; // d(eff)/dp1 err = sqrt(err0*err0 + err1*err1 - 2*rho*err0*err1); double avgrate = 0.6867; float wz_ele_eta[] = {0.588717,0.603674,0.666135,0.747493,0.762202,0.675051,0.751606,0.745569,0.665333,0.610432,0.592693,}; float ewz_ele_eta[] ={0.00292902,0.002476,0.00241209,0.00182319,0.00194339,0.00299785,0.00197339,0.00182004,0.00241793,0.00245997,0.00290394,}; int ibin = 3; if (eta >= -2.5 && eta < -2.0) ibin = 0; if (eta >= -2.0 && eta < -1.5) ibin = 1; if (eta >= -1.5 && eta < -1.0) ibin = 2; if (eta >= -1.0 && eta < -0.5) ibin = 3; if (eta >= -0.5 && eta < -0.1) ibin = 4; if (eta >= -0.1 && eta < 0.1) ibin = 5; if (eta >= 0.1 && eta < 0.5) ibin = 6; if (eta >= 0.5 && eta < 1.0) ibin = 7; if (eta >= 1.0 && eta < 1.5) ibin = 8; if (eta >= 1.5 && eta < 2.0) ibin = 9; if (eta >= 2.0 && eta < 2.5) ibin = 10; double eff_eta = wz_ele_eta[ibin]; double err_eta = ewz_ele_eta[ibin]; eff = (eff*eff_eta)/avgrate; } if (flavor == 12) { // weight electron from tau float rho = 0.884; float p0 = 6.799; float p1 = 0.842; float ep0= 0.664; float ep1= 0.016; eff = p1 - p0/pt; double err0 = ep0/pt; // d(eff)/dp0 double err1 = ep1; // d(eff)/dp1 err = sqrt(err0*err0 + err1*err1 - 2*rho*err0*err1); double avgrate = 0.5319; float wz_elet_eta[] = {0.468945,0.465953,0.489545,0.58709,0.59669,0.515829,0.59284,0.575828,0.498181,0.463536,0.481738,}; float ewz_elet_eta[] ={0.00933795,0.00780868,0.00792679,0.00642083,0.00692652,0.0101568,0.00698452,0.00643524,0.0080002,0.00776238,0.0094699,}; int ibin = 3; if (eta >= -2.5 && eta < -2.0) ibin = 0; if (eta >= -2.0 && eta < -1.5) ibin = 1; if (eta >= -1.5 && eta < -1.0) ibin = 2; if (eta >= -1.0 && eta < -0.5) ibin = 3; if (eta >= -0.5 && eta < -0.1) ibin = 4; if (eta >= -0.1 && eta < 0.1) ibin = 5; if (eta >= 0.1 && eta < 0.5) ibin = 6; if (eta >= 0.5 && eta < 1.0) ibin = 7; if (eta >= 1.0 && eta < 1.5) ibin = 8; if (eta >= 1.5 && eta < 2.0) ibin = 9; if (eta >= 2.0 && eta < 2.5) ibin = 10; double eff_eta = wz_elet_eta[ibin]; double err_eta = ewz_elet_eta[ibin]; eff = (eff*eff_eta)/avgrate; } if (flavor == 13) {// weight prompt muon //if eta>0.1 float p0 = -18.21; float p1 = 14.83; float p2 = 0.9312; float ep0= 5.06; float ep1= 1.9; float ep2=0.00069; if ( fabs(eta) < 0.1) { p0 = 7.459; p1 = 2.615; p2 = 0.5138; ep0 = 10.4; ep1 = 4.934; ep2 = 0.0034; } double arg = ( pt-p0 )/( 2.*p1 ) ; eff = 0.5 * p2 * (1.+erf(arg)); err = 0.1*eff; } if (flavor == 14) {// weight muon from tau if (fabs(eta) < 0.1) { float p0 = -1.756; float p1 = 12.38; float p2 = 0.4441; //float ep0= 10.39; float ep1= 7.9; float ep2=0.022; double arg = ( pt-p0 )/( 2.*p1 ) ; eff = 0.5 * p2 * (1.+erf(arg)); err = 0.1*eff; } else { float p0 = 2.102; float p1 = 0.8293; //float ep0= 0.271; float ep1= 0.0083; eff = p1 - p0/pt; //double err0 = ep0/pt; // d(eff)/dp0 //double err1 = ep1; // d(eff)/dp1 //err = sqrt(err0*err0 + err1*err1 - 2*rho*err0*err1); } } if (flavor == 15) {// weight hadronic tau 1p float wz_tau1p[] = {0.0249278,0.146978,0.225049,0.229212,0.21519,0.206152,0.201559,0.197917,0.209249,0.228336,0.193548,}; float ewz_tau1p[] ={0.00178577,0.00425252,0.00535052,0.00592126,0.00484684,0.00612941,0.00792099,0.0083006,0.0138307,0.015568,0.0501751,}; int ibin = 0; if (pt > 15) ibin = 1; if (pt > 20) ibin = 2; if (pt > 25) ibin = 3; if (pt > 30) ibin = 4; if (pt > 40) ibin = 5; if (pt > 50) ibin = 6; if (pt > 60) ibin = 7; if (pt > 80) ibin = 8; if (pt > 100) ibin = 9; if (pt > 200) ibin = 10; eff = wz_tau1p[ibin]; err = ewz_tau1p[ibin]; double avgrate = 0.1718; float wz_tau1p_eta[] = {0.162132,0.176393,0.139619,0.178813,0.185144,0.210027,0.203937,0.178688,0.137034,0.164216,0.163713,}; float ewz_tau1p_eta[] ={0.00706705,0.00617989,0.00506798,0.00525172,0.00581865,0.00865675,0.00599245,0.00529877,0.00506368,0.00617025,0.00726219,}; ibin = 3; if (eta >= -2.5 && eta < -2.0) ibin = 0; if (eta >= -2.0 && eta < -1.5) ibin = 1; if (eta >= -1.5 && eta < -1.0) ibin = 2; if (eta >= -1.0 && eta < -0.5) ibin = 3; if (eta >= -0.5 && eta < -0.1) ibin = 4; if (eta >= -0.1 && eta < 0.1) ibin = 5; if (eta >= 0.1 && eta < 0.5) ibin = 6; if (eta >= 0.5 && eta < 1.0) ibin = 7; if (eta >= 1.0 && eta < 1.5) ibin = 8; if (eta >= 1.5 && eta < 2.0) ibin = 9; if (eta >= 2.0 && eta < 2.5) ibin = 10; double eff_eta = wz_tau1p_eta[ibin]; double err_eta = ewz_tau1p_eta[ibin]; eff = (eff*eff_eta)/avgrate; } if (flavor == 16) { //weight hadronic tau 3p float wz_tau3p[] = {0.000587199,0.00247181,0.0013031,0.00280112,}; float ewz_tau3p[] ={0.000415091,0.000617187,0.000582385,0.00197792,}; int ibin = 0; if (pt > 15) ibin = 1; if (pt > 20) ibin = 2; if (pt > 40) ibin = 3; if (pt > 80) ibin = 4; eff = wz_tau3p[ibin]; err = ewz_tau3p[ibin]; } return eff; } /// function giving observed upper limit (visible cross-section) double Get_Upper_Limit(std::string Signal_Region, bool Observed) { std::map Upper_Limits_Observed; std::map Upper_Limits_Expected; Upper_Limits_Observed["HTlep_3l_offZ_cut_0"] = 11.; Upper_Limits_Observed["HTlep_3l_offZ_cut_100"] = 8.7; Upper_Limits_Observed["HTlep_3l_offZ_cut_150"] = 4.0; Upper_Limits_Observed["HTlep_3l_offZ_cut_200"] = 4.4; Upper_Limits_Observed["HTlep_3l_offZ_cut_300"] = 1.6; Upper_Limits_Observed["HTlep_2ltau_offZ_cut_0"] = 25.; Upper_Limits_Observed["HTlep_2ltau_offZ_cut_100"] = 14.; Upper_Limits_Observed["HTlep_2ltau_offZ_cut_150"] = 6.1; Upper_Limits_Observed["HTlep_2ltau_offZ_cut_200"] = 3.3; Upper_Limits_Observed["HTlep_2ltau_offZ_cut_300"] = 1.2; Upper_Limits_Observed["HTlep_3l_onZ_cut_0"] = 48.; Upper_Limits_Observed["HTlep_3l_onZ_cut_100"] = 38.; Upper_Limits_Observed["HTlep_3l_onZ_cut_150"] = 14.; Upper_Limits_Observed["HTlep_3l_onZ_cut_200"] = 7.2; Upper_Limits_Observed["HTlep_3l_onZ_cut_300"] = 4.5; Upper_Limits_Observed["HTlep_2ltau_onZ_cut_0"] = 85.; Upper_Limits_Observed["HTlep_2ltau_onZ_cut_100"] = 53.; Upper_Limits_Observed["HTlep_2ltau_onZ_cut_150"] = 11.0; Upper_Limits_Observed["HTlep_2ltau_onZ_cut_200"] = 5.2; Upper_Limits_Observed["HTlep_2ltau_onZ_cut_300"] = 3.0; Upper_Limits_Observed["METStrong_3l_offZ_cut_0"] = 2.6; Upper_Limits_Observed["METStrong_3l_offZ_cut_50"] = 2.1; Upper_Limits_Observed["METStrong_3l_offZ_cut_75"] = 2.1; Upper_Limits_Observed["METStrong_2ltau_offZ_cut_0"] = 4.2; Upper_Limits_Observed["METStrong_2ltau_offZ_cut_50"] = 3.1; Upper_Limits_Observed["METStrong_2ltau_offZ_cut_75"] = 2.6; Upper_Limits_Observed["METStrong_3l_onZ_cut_20"] = 11.0; Upper_Limits_Observed["METStrong_3l_onZ_cut_50"] = 6.4; Upper_Limits_Observed["METStrong_3l_onZ_cut_75"] = 5.1; Upper_Limits_Observed["METStrong_2ltau_onZ_cut_20"] = 5.9; Upper_Limits_Observed["METStrong_2ltau_onZ_cut_50"] = 3.4; Upper_Limits_Observed["METStrong_2ltau_onZ_cut_75"] = 1.2; Upper_Limits_Observed["METWeak_3l_offZ_cut_0"] = 11.; Upper_Limits_Observed["METWeak_3l_offZ_cut_50"] = 5.3; Upper_Limits_Observed["METWeak_3l_offZ_cut_75"] = 3.1; Upper_Limits_Observed["METWeak_2ltau_offZ_cut_0"] = 23.; Upper_Limits_Observed["METWeak_2ltau_offZ_cut_50"] = 4.3; Upper_Limits_Observed["METWeak_2ltau_offZ_cut_75"] = 3.1; Upper_Limits_Observed["METWeak_3l_onZ_cut_20"] = 41.; Upper_Limits_Observed["METWeak_3l_onZ_cut_50"] = 16.; Upper_Limits_Observed["METWeak_3l_onZ_cut_75"] = 8.0; Upper_Limits_Observed["METWeak_2ltau_onZ_cut_20"] = 80.; Upper_Limits_Observed["METWeak_2ltau_onZ_cut_50"] = 4.4; Upper_Limits_Observed["METWeak_2ltau_onZ_cut_75"] = 1.8; Upper_Limits_Observed["Meff_3l_offZ_cut_0"] = 11.; Upper_Limits_Observed["Meff_3l_offZ_cut_150"] = 8.1; Upper_Limits_Observed["Meff_3l_offZ_cut_300"] = 3.1; Upper_Limits_Observed["Meff_3l_offZ_cut_500"] = 2.1; Upper_Limits_Observed["Meff_2ltau_offZ_cut_0"] = 25.; Upper_Limits_Observed["Meff_2ltau_offZ_cut_150"] = 12.; Upper_Limits_Observed["Meff_2ltau_offZ_cut_300"] = 3.9; Upper_Limits_Observed["Meff_2ltau_offZ_cut_500"] = 2.2; Upper_Limits_Observed["Meff_3l_onZ_cut_0"] = 48.; Upper_Limits_Observed["Meff_3l_onZ_cut_150"] = 37.; Upper_Limits_Observed["Meff_3l_onZ_cut_300"] = 11.; Upper_Limits_Observed["Meff_3l_onZ_cut_500"] = 4.8; Upper_Limits_Observed["Meff_2ltau_onZ_cut_0"] = 85.; Upper_Limits_Observed["Meff_2ltau_onZ_cut_150"] = 28.; Upper_Limits_Observed["Meff_2ltau_onZ_cut_300"] = 5.9; Upper_Limits_Observed["Meff_2ltau_onZ_cut_500"] = 1.9; Upper_Limits_Observed["MeffStrong_3l_offZ_cut_0"] = 3.8; Upper_Limits_Observed["MeffStrong_3l_offZ_cut_150"] = 3.8; Upper_Limits_Observed["MeffStrong_3l_offZ_cut_300"] = 2.8; Upper_Limits_Observed["MeffStrong_3l_offZ_cut_500"] = 2.1; Upper_Limits_Observed["MeffStrong_2ltau_offZ_cut_0"] = 3.9; Upper_Limits_Observed["MeffStrong_2ltau_offZ_cut_150"] = 4.0; Upper_Limits_Observed["MeffStrong_2ltau_offZ_cut_300"] = 2.9; Upper_Limits_Observed["MeffStrong_2ltau_offZ_cut_500"] = 1.5; Upper_Limits_Observed["MeffStrong_3l_onZ_cut_0"] = 10.0; Upper_Limits_Observed["MeffStrong_3l_onZ_cut_150"] = 10.0; Upper_Limits_Observed["MeffStrong_3l_onZ_cut_300"] = 6.8; Upper_Limits_Observed["MeffStrong_3l_onZ_cut_500"] = 3.9; Upper_Limits_Observed["MeffStrong_2ltau_onZ_cut_0"] = 1.6; Upper_Limits_Observed["MeffStrong_2ltau_onZ_cut_150"] = 1.4; Upper_Limits_Observed["MeffStrong_2ltau_onZ_cut_300"] = 1.5; Upper_Limits_Observed["MeffStrong_2ltau_onZ_cut_500"] = 0.9; // expected upper limits are also given but not used in this Analysis Upper_Limits_Expected["HTlep_3l_offZ_cut_0"] = 11.; Upper_Limits_Expected["HTlep_3l_offZ_cut_100"] = 8.5; Upper_Limits_Expected["HTlep_3l_offZ_cut_150"] = 4.6; Upper_Limits_Expected["HTlep_3l_offZ_cut_200"] = 3.6; Upper_Limits_Expected["HTlep_3l_offZ_cut_300"] = 1.9; Upper_Limits_Expected["HTlep_2ltau_offZ_cut_0"] = 23.; Upper_Limits_Expected["HTlep_2ltau_offZ_cut_100"] = 14.; Upper_Limits_Expected["HTlep_2ltau_offZ_cut_150"] = 6.4; Upper_Limits_Expected["HTlep_2ltau_offZ_cut_200"] = 3.6; Upper_Limits_Expected["HTlep_2ltau_offZ_cut_300"] = 1.5; Upper_Limits_Expected["HTlep_3l_onZ_cut_0"] = 33.; Upper_Limits_Expected["HTlep_3l_onZ_cut_100"] = 25.; Upper_Limits_Expected["HTlep_3l_onZ_cut_150"] = 12.; Upper_Limits_Expected["HTlep_3l_onZ_cut_200"] = 6.5; Upper_Limits_Expected["HTlep_3l_onZ_cut_300"] = 3.1; Upper_Limits_Expected["HTlep_2ltau_onZ_cut_0"] = 94.; Upper_Limits_Expected["HTlep_2ltau_onZ_cut_100"] = 61.; Upper_Limits_Expected["HTlep_2ltau_onZ_cut_150"] = 9.9; Upper_Limits_Expected["HTlep_2ltau_onZ_cut_200"] = 4.5; Upper_Limits_Expected["HTlep_2ltau_onZ_cut_300"] = 1.9; Upper_Limits_Expected["METStrong_3l_offZ_cut_0"] = 3.1; Upper_Limits_Expected["METStrong_3l_offZ_cut_50"] = 2.4; Upper_Limits_Expected["METStrong_3l_offZ_cut_75"] = 2.3; Upper_Limits_Expected["METStrong_2ltau_offZ_cut_0"] = 4.8; Upper_Limits_Expected["METStrong_2ltau_offZ_cut_50"] = 3.3; Upper_Limits_Expected["METStrong_2ltau_offZ_cut_75"] = 2.1; Upper_Limits_Expected["METStrong_3l_onZ_cut_20"] = 8.7; Upper_Limits_Expected["METStrong_3l_onZ_cut_50"] = 4.9; Upper_Limits_Expected["METStrong_3l_onZ_cut_75"] = 3.8; Upper_Limits_Expected["METStrong_2ltau_onZ_cut_20"] = 7.3; Upper_Limits_Expected["METStrong_2ltau_onZ_cut_50"] = 2.8; Upper_Limits_Expected["METStrong_2ltau_onZ_cut_75"] = 1.5; Upper_Limits_Expected["METWeak_3l_offZ_cut_0"] = 10.; Upper_Limits_Expected["METWeak_3l_offZ_cut_50"] = 4.7; Upper_Limits_Expected["METWeak_3l_offZ_cut_75"] = 3.0; Upper_Limits_Expected["METWeak_2ltau_offZ_cut_0"] = 21.; Upper_Limits_Expected["METWeak_2ltau_offZ_cut_50"] = 4.0; Upper_Limits_Expected["METWeak_2ltau_offZ_cut_75"] = 2.6; Upper_Limits_Expected["METWeak_3l_onZ_cut_20"] = 30.; Upper_Limits_Expected["METWeak_3l_onZ_cut_50"] = 10.; Upper_Limits_Expected["METWeak_3l_onZ_cut_75"] = 5.4; Upper_Limits_Expected["METWeak_2ltau_onZ_cut_20"] = 88.; Upper_Limits_Expected["METWeak_2ltau_onZ_cut_50"] = 5.5; Upper_Limits_Expected["METWeak_2ltau_onZ_cut_75"] = 2.2; Upper_Limits_Expected["Meff_3l_offZ_cut_0"] = 11.; Upper_Limits_Expected["Meff_3l_offZ_cut_150"] = 8.8; Upper_Limits_Expected["Meff_3l_offZ_cut_300"] = 3.7; Upper_Limits_Expected["Meff_3l_offZ_cut_500"] = 2.1; Upper_Limits_Expected["Meff_2ltau_offZ_cut_0"] = 23.; Upper_Limits_Expected["Meff_2ltau_offZ_cut_150"] = 13.; Upper_Limits_Expected["Meff_2ltau_offZ_cut_300"] = 4.9; Upper_Limits_Expected["Meff_2ltau_offZ_cut_500"] = 2.4; Upper_Limits_Expected["Meff_3l_onZ_cut_0"] = 33.; Upper_Limits_Expected["Meff_3l_onZ_cut_150"] = 25.; Upper_Limits_Expected["Meff_3l_onZ_cut_300"] = 9.; Upper_Limits_Expected["Meff_3l_onZ_cut_500"] = 3.9; Upper_Limits_Expected["Meff_2ltau_onZ_cut_0"] = 94.; Upper_Limits_Expected["Meff_2ltau_onZ_cut_150"] = 35.; Upper_Limits_Expected["Meff_2ltau_onZ_cut_300"] = 6.8; Upper_Limits_Expected["Meff_2ltau_onZ_cut_500"] = 2.5; Upper_Limits_Expected["MeffStrong_3l_offZ_cut_0"] = 3.9; Upper_Limits_Expected["MeffStrong_3l_offZ_cut_150"] = 3.9; Upper_Limits_Expected["MeffStrong_3l_offZ_cut_300"] = 3.0; Upper_Limits_Expected["MeffStrong_3l_offZ_cut_500"] = 2.0; Upper_Limits_Expected["MeffStrong_2ltau_offZ_cut_0"] = 3.8; Upper_Limits_Expected["MeffStrong_2ltau_offZ_cut_150"] = 3.9; Upper_Limits_Expected["MeffStrong_2ltau_offZ_cut_300"] = 3.1; Upper_Limits_Expected["MeffStrong_2ltau_offZ_cut_500"] = 1.6; Upper_Limits_Expected["MeffStrong_3l_onZ_cut_0"] = 6.9; Upper_Limits_Expected["MeffStrong_3l_onZ_cut_150"] = 7.1; Upper_Limits_Expected["MeffStrong_3l_onZ_cut_300"] = 4.9; Upper_Limits_Expected["MeffStrong_3l_onZ_cut_500"] = 3.0; Upper_Limits_Expected["MeffStrong_2ltau_onZ_cut_0"] = 2.4; Upper_Limits_Expected["MeffStrong_2ltau_onZ_cut_150"] = 2.5; Upper_Limits_Expected["MeffStrong_2ltau_onZ_cut_300"] = 2.0; Upper_Limits_Expected["MeffStrong_2ltau_onZ_cut_500"] = 1.1; if (Observed) return Upper_Limits_Observed[Signal_Region]; else return Upper_Limits_Expected[Signal_Region]; } /// function checking if there is an OSSF lepton pair or a combination of 3 leptons with an invariant mass close to the Z mass int isonZ ( const ParticleVector & particles) { int onZ = 0; double best_mass_2 = 999.; double best_mass_3 = 999.; // loop over all 2 particle combinations to find invariant mass of OSSF pair closest to Z mass foreach ( const Particle & p1 , particles ) { foreach ( const Particle & p2 , particles ) { double mass_difference_2_old = fabs(91.0 - best_mass_2); double mass_difference_2_new = fabs(91.0 - (p1.momentum() + p2.momentum()).mass()/GeV); // if particle combination is OSSF pair calculate mass difference to Z mass if ( (p1.pdgId()*p2.pdgId() == -121 || p1.pdgId()*p2.pdgId() == -169) ) { // get invariant mass closest to Z mass if (mass_difference_2_new < mass_difference_2_old) best_mass_2 = (p1.momentum() + p2.momentum()).mass()/GeV; // in case there is an OSSF pair take also 3rd lepton into account (e.g. from FSR and photon to electron conversion) foreach ( const Particle & p3 , particles ) { double mass_difference_3_old = fabs(91.0 - best_mass_3); double mass_difference_3_new = fabs(91.0 - (p1.momentum() + p2.momentum() + p3.momentum()).mass()/GeV); if (mass_difference_3_new < mass_difference_3_old) best_mass_3 = (p1.momentum() + p2.momentum() + p3.momentum()).mass()/GeV; } } } } // pick the minimum invariant mass of the best OSSF pair combination and the best 3 lepton combination double best_mass = min(best_mass_2, best_mass_3); // if this mass is in a 20 GeV window around the Z mass, the event is classified as onZ if ( fabs(91.0 - best_mass) < 20. ) onZ = 1; return onZ; } }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2012_I1204447); }