// -*- 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_EXOTICS_MULTILEPTON_7TEV : public Analysis { public: /// Constructor ATLAS_EXOTICS_MULTILEPTON_7TEV() : Analysis("ATLAS_EXOTICS_MULTILEPTON_7TEV") { } /// public variables // some MC event generators use different units than others, so a factor 1000 has to be introduced sometimes, please check the histograms and adjust if necessary double Rivet2Athena; double MyMeV; double MyGeV; // counting number of total events int number_events_total; // list of signal regions and event counts per signal region std::vector Signal_Regions; std::map EventCountsPerSR; public: /// Book histograms and initialise projections before the run void init() { // some MC event generators use different units than others, so a factor 1000 has to be introduced sometimes, please check the histograms and adjust if necessary Rivet2Athena = 1.; // Rivet2Athena = 1000.; MyMeV = MeV/Rivet2Athena; MyGeV = GeV/Rivet2Athena; // counting number of total events number_events_total =0; // 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; /// final state including all charged and neutral particles const FinalState fs(-5.0, 5.0, 1000*MyMeV); addProjection(fs, "FS"); /// final state including all charged particles addProjection(ChargedFinalState(-2.5, 2.5, 1000*MyMeV), "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*MyGeV),"UFS"); /// final state including all electrons IdentifiedFinalState elecs(-2.47,2.47, 10.0*MyGeV); elecs.acceptIdPair(ELECTRON); addProjection(elecs, "elecs"); /// final state including all muons IdentifiedFinalState muons(-2.5,2.5, 10.0*MyGeV); 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(); number_events_total +=1; /// 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 = apply_reco_eff(muon_id,mu.momentum().pT()/MyGeV,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 = apply_reco_eff(elec_id,e.momentum().pT()/MyGeV,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()/MyGeV <= 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 = apply_reco_eff(tau_id,tau.momentum().pT()/MyGeV,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*MyGeV) ) { 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()/MyGeV; /// overlap removal /// electron - electron ParticleVector electron_candidates_2; for(unsigned int ie=0;ie(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()/MyGeV>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()/MyGeV>25.) && !( !trigger_mu.empty() && trigger_mu[0].momentum().perp()/MyGeV>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()/MyGeV,weight); _h_pt_2_3l->fill(recon_leptons[1].momentum().perp()/MyGeV,weight); _h_pt_3_3l->fill(recon_leptons[2].momentum().perp()/MyGeV,weight); HTlep = recon_leptons[0].momentum().perp()/MyGeV + recon_leptons[1].momentum().perp()/MyGeV + recon_leptons[2].momentum().perp()/MyGeV; 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()/MyGeV,weight); _h_pt_2_2ltau->fill(recon_leptons[1].momentum().perp()/MyGeV,weight); _h_pt_3_2ltau->fill(recon_tau[0].momentum().perp()/MyGeV,weight); HTlep = recon_leptons[0].momentum().perp()/MyGeV + recon_leptons[1].momentum().perp()/MyGeV + recon_tau[0].momentum().perp()/MyGeV; 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()/MyGeV; /// calculate meff double meff = eTmiss + HTjets; ParticleVector all_leptons; foreach ( const Particle & e , recon_e ) { meff += e.momentum().perp()/MyGeV; all_leptons.push_back( e ); } foreach ( const Particle & mu, recon_mu ) { meff += mu.momentum().perp()/MyGeV; all_leptons.push_back( mu ); } foreach ( const Particle & tau, recon_tau ) { meff += tau.momentum().perp()/MyGeV; 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); } /// Normalise histograms etc., after the run void finalize() { cout << "Number of total events: " << number_events_total << endl; // normalize to an integrated luminosity of 1 fb-1 double norm = crossSection()/femtobarn/sumOfWeights(); // check if cross-section is 0, in this case the right cross-section should be given by hand (or the inputfile/format should be adjusted) if (norm == 0) { cout << "-----------------------------------------------------------------------------------------------------------" << endl; cout << "ERROR: cross-section is 0! Therefore the number of events is normalized to 0! \n Either the input file is broken or an old HepMC version is used. The right cross-section can be put in by hand here." << endl; cout << "-----------------------------------------------------------------------------------------------------------" << endl; double right_cross_section = 0.; // please enter right cross section in pb norm = right_cross_section*1000./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 << "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; if (signal_events_best_SR > exp_UL_best_SR) { cout << "Using the EXPECTED limits (visible cross-section) of the analysis: " << endl; 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; if (signal_events_best_SR > obs_UL_best_SR) { cout << "Using the OBSERVED limits (visible cross-section) of the analysis: " << endl; 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, crossSection()/femtobarn/sumOfWeights()); scale(_h_HTjets_all, crossSection()/femtobarn/sumOfWeights()); scale(_h_MET_all, crossSection()/femtobarn/sumOfWeights()); scale(_h_Meff_all, crossSection()/femtobarn/sumOfWeights()); scale(_h_pt_1_3l, crossSection()/femtobarn/sumOfWeights()); scale(_h_pt_2_3l, crossSection()/femtobarn/sumOfWeights()); scale(_h_pt_3_3l, crossSection()/femtobarn/sumOfWeights()); scale(_h_pt_1_2ltau, crossSection()/femtobarn/sumOfWeights()); scale(_h_pt_2_2ltau, crossSection()/femtobarn/sumOfWeights()); scale(_h_pt_3_2ltau, crossSection()/femtobarn/sumOfWeights()); scale(_h_e_n, crossSection()/femtobarn/sumOfWeights()); scale(_h_mu_n, crossSection()/femtobarn/sumOfWeights()); scale(_h_tau_n, crossSection()/femtobarn/sumOfWeights()); scale(_h_excluded, crossSection()/femtobarn/sumOfWeights()); } } 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 cut_values = Get_Cuts_Per_Signal_Region(kinematic_variables.at(i0),onZ); // loop over all cut values for(unsigned int i2=0; i2 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); // the MET < 20 GeV onZ region is a control region 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) { // 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.at(i)) EventCountsPerSR[("HTlep_" + basic_signal_region + "_cut_" + int2str(cut_values.at(i))).c_str()]++; } // 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.at(i) && HTjets > 100.) EventCountsPerSR[("METStrong_" + basic_signal_region + "_cut_" + int2str(cut_values.at(i))).c_str()]++; } // 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.at(i) && HTjets <= 100.) EventCountsPerSR[("METWeak_" + basic_signal_region + "_cut_" + int2str(cut_values.at(i))).c_str()]++; } // 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.at(i)) EventCountsPerSR[("Meff_" + basic_signal_region + "_cut_" + int2str(cut_values.at(i))).c_str()]++; } // 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.at(i) && eTmiss > 75.) EventCountsPerSR[("MeffStrong_" + basic_signal_region + "_cut_" + int2str(cut_values.at(i))).c_str()]++; } } /// 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 = 0.; double err_eta = 0; eff_eta = wz_ele_eta[ibin]; 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 = 0.; double err_eta = 0; eff_eta = wz_elet_eta[ibin]; 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 = 0.; double err_eta = 0; eff_eta = wz_tau1p_eta[ibin]; 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()/MyGeV); // 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()/MyGeV; // 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()/MyGeV); if (mass_difference_3_new < mass_difference_3_old) best_mass_3 = (p1.momentum()+p2.momentum()+p3.momentum()).mass()/MyGeV; } } } } // 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_EXOTICS_MULTILEPTON_7TEV); }