|
[Rivet-svn] r3275 - in trunk: data/anainfo src/Analysesblackhole at projects.hepforge.org blackhole at projects.hepforge.orgTue Aug 2 16:27:14 BST 2011
Author: dgrell Date: Tue Aug 2 16:27:14 2011 New Revision: 3275 Log: Two more SUSY analyses from Angela Chen. Unvalidated. Added: trunk/data/anainfo/ATLAS_2011_CONF_2011_098.info trunk/data/anainfo/ATLAS_2011_S9041966.info trunk/src/Analyses/ATLAS_2011_CONF_2011_098.cc trunk/src/Analyses/ATLAS_2011_S9041966.cc Modified: trunk/data/anainfo/Makefile.am trunk/src/Analyses/Makefile.am Added: trunk/data/anainfo/ATLAS_2011_CONF_2011_098.info ============================================================================== --- /dev/null 00:00:00 1970 (empty, because file is newly added) +++ trunk/data/anainfo/ATLAS_2011_CONF_2011_098.info Tue Aug 2 16:27:14 2011 (r3275) @@ -0,0 +1,24 @@ +Name: ATLAS_2011_CONF_2011_098 +Year: 2011 +Summary: B-jets Search for Supersymmetry with 0-leptons +Experiment: ATLAS +Collider: LHC +SpiresID: +Status: UNVALIDATED +Authors: + - Angela Chen <aqchen at fas.harvard.edu> +References: + - arXiv:nnnn.nnnn +RunInfo: + BSM signal events at 7000 GeV. +NumEvents: 25000 for BSM signals +Beams: [p+, p+] +Energies: [7000] +Description: + 'Search for supersymmmetric particles by ATLAS at 7 TeV in events with b-jets, large missing + energy, and no leptons. + Event counts in four signal regions (1 b-jet, m_eff>500*GeV; 1 b-jet, m_eff>700*GeV; 2 b-jets, + m_eff>500*GeV; 2 b-jets, m_eff>700*GeV) are implemented as one-bin histograms. + Histograms for missing transverse energy, effective mass, and pT of the leading jet are + implemented for the 1 b-tag and 2 b-tag signal regions.' + Added: trunk/data/anainfo/ATLAS_2011_S9041966.info ============================================================================== --- /dev/null 00:00:00 1970 (empty, because file is newly added) +++ trunk/data/anainfo/ATLAS_2011_S9041966.info Tue Aug 2 16:27:14 2011 (r3275) @@ -0,0 +1,22 @@ +Name: ATLAS_2011_S9041966 +Year: 2011 +Summary: 1-lepton and 2-lepton search for first or second generation leptoquarks +Experiment: ATLAS +Collider: LHC +SpiresID: 9041966 +Status: UNVALIDATED +Authors: + - Angela Chen <aqchen at fas.harvard.edu> +References: + - arXiv:1104.4481 +RunInfo: + BSM signal events at 7000 GeV. +NumEvents: 25000 for BSM signals +Beams: [p+, p+] +Energies: [7000] +Description: + 'Single and dilepton search for first and second generation scalar leptoquarks by ATLAS at 7 TeV. + Event counts in four signal regions (single lepton and dilepton for first and second generation) are implemented as one-bin histograms. + Histograms for event transverse energy are implemented for dilepton signal regions and histograms for leptoquark mass are implemented for single lepton signal regions. + Histograms for observables in six control regions are implemented.' + Modified: trunk/data/anainfo/Makefile.am ============================================================================== --- trunk/data/anainfo/Makefile.am Tue Aug 2 13:53:41 2011 (r3274) +++ trunk/data/anainfo/Makefile.am Tue Aug 2 16:27:14 2011 (r3275) @@ -19,7 +19,9 @@ ATLAS_2011_S8994773.info \ ATLAS_2011_S9002537.info \ ATLAS_2011_S9019561.info \ + ATLAS_2011_S9041966.info \ ATLAS_2011_CONF_2011_090.info \ + ATLAS_2011_CONF_2011_098.info \ ATLAS_2011_S9120807.info \ BELLE_2006_S6265367.info \ CDF_1988_S1865951.info \ Added: trunk/src/Analyses/ATLAS_2011_CONF_2011_098.cc ============================================================================== --- /dev/null 00:00:00 1970 (empty, because file is newly added) +++ trunk/src/Analyses/ATLAS_2011_CONF_2011_098.cc Tue Aug 2 16:27:14 2011 (r3275) @@ -0,0 +1,344 @@ +// -*- C++ -*- +#include "Rivet/Analysis.hh" +#include "Rivet/Tools/BinnedHistogram.hh" +#include "Rivet/RivetAIDA.hh" +#include "Rivet/Tools/Logging.hh" +#include "Rivet/Projections/FinalState.hh" +#include "Rivet/Projections/ChargedFinalState.hh" +#include "Rivet/Projections/VisibleFinalState.hh" +#include "Rivet/Projections/IdentifiedFinalState.hh" +#include "Rivet/Projections/FastJets.hh" +// #include "Rivet/Tools/RivetMT2.hh" + +namespace Rivet { + + + class ATLAS_2011_CONF_2011_098 : public Analysis { + public: + + /// @name Constructors etc. + //@{ + + /// Constructor + ATLAS_2011_CONF_2011_098() + : Analysis("ATLAS_2011_CONF_2011_098"), + +//debug variables +threeJA(0), threeJB(0), threeJC(0), threeJD(0), bj(0), jets(0), zerolept(0), eTmisscut(0) + + + { + /// Set whether your finalize method needs the generator cross section + setNeedsCrossSection(false); + } + + //@} + + + public: + + /// @name Analysis methods + //@{ + + /// Book histograms and initialise projections before the run + void init() { + + // projection to find the electrons + std::vector<std::pair<double, double> > eta_e; + eta_e.push_back(make_pair(-2.47,2.47)); + IdentifiedFinalState elecs(eta_e, 20.0*GeV); + elecs.acceptIdPair(ELECTRON); + addProjection(elecs, "elecs"); + + + // projection to find the muons + std::vector<std::pair<double, double> > eta_m; + eta_m.push_back(make_pair(-2.4,2.4)); + IdentifiedFinalState muons(eta_m, 10.0*GeV); + muons.acceptIdPair(MUON); + addProjection(muons, "muons"); + + + + /// Jet finder + addProjection(FastJets(FinalState(), FastJets::ANTIKT, 0.4), + "AntiKtJets04"); + + + // all tracks (to do deltaR with leptons) + addProjection(ChargedFinalState(-3.0,3.0),"cfs"); + + // for pTmiss + addProjection(VisibleFinalState(-4.9,4.9),"vfs"); + + + /// Book histograms + _count_threeJA = bookHistogram1D("count_threeJA", 1, 0., 1.); + _count_threeJB = bookHistogram1D("count_threeJB", 1, 0., 1.); + _count_threeJC = bookHistogram1D("count_threeJC", 1, 0., 1.); + _count_threeJD = bookHistogram1D("count_threeJD", 1, 0., 1.); + _hist_meff_1bjet = bookHistogram1D("meff_1bjet", 32, 0., 1600.); + _hist_eTmiss_1bjet = bookHistogram1D("eTmiss_1bjet", 6, 0., 600.); + _hist_pTj_1bjet = bookHistogram1D("pTjet_1bjet", 20, 0., 800.); + _hist_meff_2bjet = bookHistogram1D("meff_2bjet", 32, 0., 1600.); + _hist_eTmiss_2bjet = bookHistogram1D("eTmiss_2bjet", 6, 0., 600.); + _hist_pTj_2bjet = bookHistogram1D("pTjet_2bjet", 20, 0., 800.); + + + } + + + /// Perform the per-event analysis + void analyze(const Event& event) { + + const double weight = event.weight(); + + // Temp: calorimeter module failure with 10% acceptance loss; + // region unknown ==> randomly choose 10% of events to be vetoed + if ( rand()/RAND_MAX < 0.1 ) + vetoEvent; + + Jets tmp_cand_jets; + foreach (const Jet& jet, + applyProjection<FastJets>(event, "AntiKtJets04").jetsByPt(20.0*GeV) ) { + if ( fabs( jet.momentum().eta() ) < 2.8 ) { + tmp_cand_jets.push_back(jet); + } + } + + ParticleVector cand_e = + applyProjection<IdentifiedFinalState>(event, "elecs").particlesByPt(); + ParticleVector cand_mu = + applyProjection<IdentifiedFinalState>(event, "muons").particlesByPt(); + ParticleVector chg_tracks = + applyProjection<ChargedFinalState>(event, "cfs").particles(); + +//cerr << "cand_e.size(): " << cand_e.size() << " cand_mu.size(): " << cand_mu.size() << '\n'; + + + Jets cand_jets; + foreach ( const Jet& jet, tmp_cand_jets ) { + if ( fabs( jet.momentum().eta() ) >= 2.8 ) + cand_jets.push_back( jet ); + else { + bool away_from_e = true; + foreach ( const Particle & e, cand_e ) { + if ( deltaR(e.momentum(),jet.momentum()) <= 0.2 ) { + away_from_e = false; + break; + } + } + if ( away_from_e ) + cand_jets.push_back( jet ); + } + } + + ParticleVector cand_lept; + + bool isolated_e; + foreach ( const Particle & e, cand_e ) { + isolated_e = true; + foreach ( const Jet& jet, cand_jets ) { + if ( deltaR(e.momentum(),jet.momentum()) < 0.4 ) + isolated_e = false; + } + if ( isolated_e == true ) + cand_lept.push_back( e ); + } + + + bool isolated_mu; + foreach ( const Particle & mu, cand_mu ) { + isolated_mu = true; + foreach ( const Jet& jet, cand_jets ) { + if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 ) + isolated_mu = false; + } + if ( isolated_mu == true) + cand_lept.push_back( mu ); + } + + + // pTmiss + ParticleVector vfs_particles + = applyProjection<VisibleFinalState>(event, "vfs").particles(); + FourMomentum pTmiss; + foreach ( const Particle & p, vfs_particles ) { + pTmiss -= p.momentum(); + } + double eTmiss = pTmiss.pT(); + + + // bjets + Jets bjets; + foreach (const Jet& j, cand_jets) { + if ( j.momentum().pT() > 20*GeV ) { + if (j.containsBottom()) bjets += j; + } + } + + if (bjets.empty()) { + MSG_DEBUG("No b-jet axes in acceptance"); + vetoEvent; + } + +++bj; + + + + // Jets event selection + if ( cand_jets.size() < 3 ) + vetoEvent; + if ( cand_jets[0].momentum().pT() <= 130*GeV ) + vetoEvent; + if ( cand_jets[1].momentum().pT() <= 50*GeV || + cand_jets[2].momentum().pT() <= 50*GeV ) + vetoEvent; +++jets; + + // eTmiss cut + if ( eTmiss <= 130*GeV ) + vetoEvent; + +++eTmisscut; + + // 0-lepton requirement + if ( !cand_lept.empty() ) + vetoEvent; +++zerolept; + + // m_eff cut + double m_eff = eTmiss + + cand_jets[0].momentum().pT() + + cand_jets[1].momentum().pT() + + cand_jets[2].momentum().pT(); + + if ( eTmiss / m_eff <= 0.25 ) + vetoEvent; + + + // min_dPhi + double min_dPhi = 999.999; + for ( int i = 0; i < 3; ++i ) { + double dPhi = deltaPhi( pTmiss.phi(), cand_jets[i].momentum().phi() ); + min_dPhi = min( min_dPhi, dPhi ); + } + + if ( min_dPhi <= 0.4 ) + vetoEvent; + + + + // ==================== FILL ==================== + + + // 1 bjet + if ( bjets.size() >= 1 ) { + + _hist_meff_1bjet->fill(m_eff, weight); + _hist_eTmiss_1bjet->fill(eTmiss, weight); + _hist_pTj_1bjet->fill(cand_jets[0].momentum().pT(), weight); + + // 3JA region + if ( m_eff > 200*GeV ) { +++threeJA; + _count_threeJA->fill(0.5, weight); + } + + // 3JB region + if ( m_eff > 700*GeV ) { +++threeJB; + _count_threeJB->fill(0.5, weight); + } + } + + // 2 bjets + if ( bjets.size() >= 2 ) { + + _hist_meff_2bjet->fill(m_eff, weight); + _hist_eTmiss_2bjet->fill(eTmiss, weight); + _hist_pTj_2bjet->fill(cand_jets[0].momentum().pT(), weight); + + // 3JC region + if ( m_eff > 500*GeV ) { +++threeJC; + _count_threeJC->fill(0.5, weight); + } + + // 3JD region + if ( m_eff > 700*GeV ) { +++threeJD; + _count_threeJD->fill(0.5, weight); + } + } + + + + + } + + //@} + + void finalize() { + + + scale( _hist_meff_1bjet, 50. * 830. * crossSection()/sumOfWeights() ); + scale( _hist_eTmiss_1bjet, 100. * 830. * crossSection()/sumOfWeights() ); + scale( _hist_pTj_1bjet, 40. * 830. * crossSection()/sumOfWeights() ); + scale( _hist_meff_2bjet, 50. * 830. * crossSection()/sumOfWeights() ); + scale( _hist_eTmiss_2bjet, 100. * 830. * crossSection()/sumOfWeights() ); + scale( _hist_pTj_2bjet, 40. * 830. * crossSection()/sumOfWeights() ); + + +cerr<< '\n'<<'\n' +<< "Saw " +<< bj << " events aft bjets cut, " +<< jets << " events aft jet cuts, " +<< eTmisscut << " events aft eTmiss cut, " +<< zerolept << " events after 0-lept cut. " +<< '\n' +<< threeJA << " 3JA events, " +<< threeJB << " 3JB events, " +<< threeJC << " 3JC events, " +<< threeJD << " 3JD events. " +<< '\n' +; + + } + + private: + + /// @name Histograms + //@{ + AIDA::IHistogram1D* _count_threeJA; + AIDA::IHistogram1D* _count_threeJB; + AIDA::IHistogram1D* _count_threeJC; + AIDA::IHistogram1D* _count_threeJD; + AIDA::IHistogram1D* _hist_meff_1bjet; + AIDA::IHistogram1D* _hist_eTmiss_1bjet; + AIDA::IHistogram1D* _hist_pTj_1bjet; + AIDA::IHistogram1D* _hist_meff_2bjet; + AIDA::IHistogram1D* _hist_eTmiss_2bjet; + AIDA::IHistogram1D* _hist_pTj_2bjet; + + //@} + +// debug variables +int threeJA; +int threeJB; +int threeJC; +int threeJD; +int bj; +int jets; +int zerolept; +int eTmisscut; + + }; + + + + // This global object acts as a hook for the plugin system + AnalysisBuilder<ATLAS_2011_CONF_2011_098> plugin_ATLAS_2011_CONF_2011_098; + + +} Added: trunk/src/Analyses/ATLAS_2011_S9041966.cc ============================================================================== --- /dev/null 00:00:00 1970 (empty, because file is newly added) +++ trunk/src/Analyses/ATLAS_2011_S9041966.cc Tue Aug 2 16:27:14 2011 (r3275) @@ -0,0 +1,689 @@ +// -*- C++ -*- +#include "Rivet/Analysis.hh" +#include "Rivet/Tools/BinnedHistogram.hh" +#include "Rivet/RivetAIDA.hh" +#include "Rivet/Tools/Logging.hh" +#include "Rivet/Projections/FinalState.hh" +#include "Rivet/Projections/ChargedFinalState.hh" +#include "Rivet/Projections/VisibleFinalState.hh" +#include "Rivet/Projections/IdentifiedFinalState.hh" +#include "Rivet/Projections/FastJets.hh" + +namespace Rivet { + + + class ATLAS_2011_S9041966 : public Analysis { + public: + + /// @name Constructors etc. + //@{ + + /// Constructor + + ATLAS_2011_S9041966() + : Analysis("ATLAS_2011_S9041966"), + +//DEBUG +count(0), vetoe(0), Njetscut(0), dilept(0), candmumujj(0), candeejj(0), onelept(0), eTmisscut(0), candmvjj(0), candevjj(0), mumujj(0), eejj(0), mTonelept(0), MLQonelept(0), MtLQonelept(0), Stvonelept(0), mTev(0), MLQev(0), MtLQev(0), Stvev(0), muvjj(0), evjj(0), emuvjj(0), cande(0), candmu(0), tmpmu(0), tmpe(0), mumuZCR(0), eeZCR(0), munuW2CR(0), munuttCR(0), enuW2CR(0), enuttCR(0) + + { + /// Set whether your finalize method needs the generator cross section + setNeedsCrossSection(false); + } + + //@} + + + public: + + /// @name Analysis methods + //@{ + + /// Book histograms and initialize projections before the run + void init() { + + // projection to find the electrons + std::vector<std::pair<double, double> > eta_e; + eta_e.push_back(make_pair(-2.47,2.47)); + IdentifiedFinalState elecs(eta_e, 20.0*GeV); + elecs.acceptIdPair(ELECTRON); + addProjection(elecs, "elecs"); + + + // veto region electrons + std::vector<std::pair<double, double> > eta_v_e; + eta_v_e.push_back(make_pair(-1.52,-1.35)); + eta_v_e.push_back(make_pair( 1.35, 1.52)); + IdentifiedFinalState veto_elecs(eta_v_e, 10.0*GeV); + veto_elecs.acceptIdPair(ELECTRON); + addProjection(veto_elecs, "veto_elecs"); + +///DEBUG + // projection to find all leptons + IdentifiedFinalState all_mu_e; + all_mu_e.acceptIdPair(MUON); + all_mu_e.acceptIdPair(ELECTRON); + addProjection(all_mu_e, "all_mu_e"); //debug + + + + // projection to find the muons + std::vector<std::pair<double, double> > eta_m; + eta_m.push_back(make_pair(-2.4,2.4)); + IdentifiedFinalState muons(eta_m, 20.0*GeV); + muons.acceptIdPair(MUON); + addProjection(muons, "muons"); + + + // Jet finder + VetoedFinalState vfs; + vfs.addVetoPairDetail(MUON,20*GeV,7000*GeV); + vfs.addVetoPairDetail(ELECTRON,20*GeV,7000*GeV); + addProjection(FastJets(vfs, FastJets::ANTIKT, 0.4), + "AntiKtJets04"); + + + // all tracks (to do deltaR with leptons) + addProjection(ChargedFinalState(-3.0,3.0,0.5*GeV),"cfs"); + + + // for pTmiss + addProjection(VisibleFinalState(-4.9,4.9),"vfs"); + + + /// Book histograms + _count_mumujj = bookHistogram1D("count_2muons_dijet", 1, 0., 1.); + _count_eejj = bookHistogram1D("count_2elecs_dijet", 1, 0., 1.); + _count_muvjj = bookHistogram1D("count_muon_neutrino_dijet", 1, 0., 1.); + _count_evjj = bookHistogram1D("count_elec_neutrino_dijet", 1, 0., 1.); + + _hist_St_mumu = bookHistogram1D("hist_mumujj_St", 10, 450., 1650.); + _hist_St_ee = bookHistogram1D("hist_eejj_St", 10, 450., 1650.); + _hist_MLQ_muv = bookHistogram1D("hist_munujj_MLQ", 9, 150., 600.); + _hist_MLQ_ev = bookHistogram1D("hist_enujj_MLQ", 9, 150., 600.); + + + + _hist_St_mumu_ZCR = bookHistogram1D("CR_Zjets_St_mumu", 40, 0., 800.); + _hist_St_ee_ZCR = bookHistogram1D("CR_Zjets_Stee", 40, 0., 800.); + _hist_MLQ_munu_W2CR = bookHistogram1D("CR_W2jets_MLQ_munu", 20, 0., 400.); + _hist_MLQ_enu_W2CR = bookHistogram1D("CR_W2jets_MLQ_enu", 20, 0., 400.); + _hist_MLQ_munu_ttCR = bookHistogram1D("CR_tt_MLQ_munu", 35, 0., 700.); + _hist_MLQ_enu_ttCR = bookHistogram1D("CR_tt_MLQ_enu", 35, 0., 700.); + + } + + + + /// Perform the per-event analysis + void analyze(const Event& event) { + + const double weight = event.weight(); + +///DEBUG + count +=1; cerr<< "Event " << count << '\n'; + // debug + + + ParticleVector veto_e + = applyProjection<IdentifiedFinalState>(event, "veto_elecs").particles(); + if ( ! veto_e.empty() ) { + MSG_DEBUG("electrons in veto region"); + vetoEvent; + } +++vetoe; + + + Jets cand_jets; + foreach ( const Jet& jet, + applyProjection<FastJets>(event, "AntiKtJets04").jetsByPt(20.0*GeV) ) { + if ( fabs( jet.momentum().eta() ) < 2.8 ) { + cand_jets.push_back(jet); + } + } + + + ParticleVector candtemp_e = + applyProjection<IdentifiedFinalState>(event, "elecs").particlesByPt(); + ParticleVector candtemp_mu = + applyProjection<IdentifiedFinalState>(event,"muons").particlesByPt(); + ParticleVector cand_mu; + ParticleVector cand_e; + ParticleVector vfs_particles + = applyProjection<VisibleFinalState>(event, "vfs").particles(); + + + // pTcone around muon track + foreach ( const Particle & mu, candtemp_mu ) { +++tmpmu; + double pTinCone = -mu.momentum().pT(); + foreach ( const Particle & track, vfs_particles ) { + if ( deltaR(mu.momentum(),track.momentum()) < 0.2 ) + pTinCone += track.momentum().pT(); + } + if ( pTinCone/mu.momentum().pT() < 0.25 ) +++candmu; + cand_mu.push_back(mu); + } + + // pTcone around electron + foreach ( const Particle e, candtemp_e ) { +++tmpe; + double pTinCone = -e.momentum().pT(); + foreach ( const Particle & track, vfs_particles ) { + if ( deltaR(e.momentum(),track.momentum()) < 0.2 ) + pTinCone += track.momentum().pT(); + } + if ( pTinCone/e.momentum().pT() < 0.2 ) +++cande; + cand_e.push_back(e); + } + + if ( cand_e.empty() && cand_mu.empty() ) { +cerr<<" ->Event vetoed. No candidate lept"<<'\n'; + vetoEvent; + } + + +//DEBUG +else{ +foreach (const Particle & mu, cand_mu) { +cerr << "cand mu: " << "Id " << mu.pdgId() << " eta " << mu.momentum().eta() << " pT " << mu.momentum().pT() << '\n'; +} +foreach (const Particle & lepton, cand_e) { +cerr << "cand e: " << "Id " << lepton.pdgId() << " eta " << lepton.momentum().eta() << " pT " << lepton.momentum().pT() << '\n'; +}} // debug + + + + // pTmiss + FourMomentum pTmiss; + foreach ( const Particle & p, vfs_particles ) { + pTmiss -= p.momentum(); + } + double eTmiss = pTmiss.pT(); + + + // discard jets that overlap with leptons + Jets recon_jets; + foreach ( const Jet& jet, cand_jets ) { + bool away_from_lept = true; + foreach ( const Particle e, cand_e ) { + if ( deltaR(e.momentum(),jet.momentum()) <= 0.5 ) { + away_from_lept = false; + break; + } + } + foreach ( const Particle & mu, cand_mu ) { + if ( deltaR(mu.momentum(),jet.momentum()) <= 0.5 ) { + away_from_lept = false; + break; + } + } + if ( away_from_lept ) + recon_jets.push_back( jet ); + } + + + +//DEBUG +cerr << " Num of recon jets: " << recon_jets.size() << '\n'; +cerr << " Num of cand e: " << cand_e.size() << '\n'; +cerr << " Num of cand mu: " << cand_mu.size() << '\n'; +//debug + + + + // ================ OBSERVABLES ================ + + + // At least 2 hard jets + if ( recon_jets.size() < 2 ) { +cerr << " ->Event vetoed. Not enough hard jets." << '\n'; + vetoEvent; + } +++Njetscut; + + + // Initialize variables for observables + double M_ll, M_LQ, St_ll, Mt_LQ, St_v, mT; + FourMomentum p_l, p_l1, p_l2, p_j[2]; + p_j[0] = recon_jets[0].momentum(); + p_j[1] = recon_jets[1].momentum(); + + ParticleVector dilept_pair; + bool single_lept = false; + + if ( cand_mu.size() == 2 && cand_e.empty() ) { +++candmumujj; + foreach ( const Particle& mu, cand_mu ) + dilept_pair.push_back(mu); + } + else if ( cand_e.size() == 2 && cand_mu.empty() ) { +++candeejj; + foreach ( const Particle& e, cand_e ) + dilept_pair.push_back(e); + } + else if ( cand_mu.size() == 1 && cand_e.empty() ) { +++candmvjj; + p_l = cand_mu[0].momentum(); + single_lept = true; + } + else if ( cand_e.size() == 1 && cand_mu.empty() ) { +++candevjj; + p_l = cand_e[0].momentum(); + single_lept = true; + } + + // Dilepton channel observables + if ( ! dilept_pair.empty() ) { + + double E_l1, E_l2, E_j1, E_j2; + double tmpM_LQ1[2], tmpM_LQ2[2], M_LQDiff1, M_LQDiff2; + + p_l1 = dilept_pair[0].momentum(); + p_l2 = dilept_pair[1].momentum(); + E_l1 = p_l1.E(); + E_l2 = p_l2.E(); + + E_j1 = p_j[0].E(); + E_j2 = p_j[1].E(); + + // Calculate possible leptoquark mass M_LQ and reconstruct average M_LQ + + tmpM_LQ1[0] = E_l1 + E_j1; + tmpM_LQ1[1] = E_l2 + E_j2; + M_LQDiff1 = abs( tmpM_LQ1[0] - tmpM_LQ1[1] ); + + tmpM_LQ2[0] = E_l1 + E_j2; + tmpM_LQ2[1] = E_l2 + E_j1; + M_LQDiff2 = abs( tmpM_LQ2[0] - tmpM_LQ2[1] ); + + if ( M_LQDiff1 > M_LQDiff2 ) + M_LQ = ( tmpM_LQ2[0] + tmpM_LQ2[1] ) / 2; + else + M_LQ = ( tmpM_LQ1[0] + tmpM_LQ1[1] ) / 2; + + // Calculate event transverse energy St + St_ll = p_l1.pT() + p_l2.pT() + p_j[0].pT() + p_j[1].pT(); + + // Dilept pair invariant mass M_ll + M_ll = E_l1 + E_l2; + + } + + // 1-lepton channel observables + else if ( single_lept ) { + + double tmpM_LQ[2], tmpMt_LQ[2], dPhi_j[2], M_LQDiff1, M_LQDiff2; + + // List of possible M_LQ, Mt_LQ pairings + + for ( int i = 0; i < 2; ++i ) { + tmpM_LQ[i] = p_l.E() + p_j[i].E(); + dPhi_j[1-i] = deltaPhi( p_j[1-i].phi(), pTmiss.phi() ); + tmpMt_LQ[i] = sqrt( 2 * p_j[1-i].pT() * eTmiss * (1 - cos( dPhi_j[1-i] )) ); + } + + // Choose pairing that gives smallest absolute difference + + M_LQDiff1 = abs( tmpM_LQ[0] - tmpMt_LQ[0] ); + M_LQDiff2 = abs( tmpM_LQ[1] - tmpMt_LQ[1] ); + + if ( M_LQDiff1 > M_LQDiff2 ) { + M_LQ = tmpM_LQ[1]; + Mt_LQ = tmpMt_LQ[1]; + } + else { + M_LQ = tmpM_LQ[0]; + Mt_LQ = tmpMt_LQ[0]; + } + + // Event transverse energy + St_v = p_l.pT() + eTmiss + p_j[0].pT() + p_j[1].pT(); + + // Transverse mass mT + double dPhi_l = deltaPhi( p_l.phi(), pTmiss.phi()); + mT = sqrt( 2 * p_l.pT() * eTmiss * (1 - cos(dPhi_l)) ); + + } + + + // ============== CONTROL REGIONS =============== + + // mumujj, Z control region + if ( cand_mu.size() == 2 ) { + if ( M_ll >= 81*GeV && M_ll <= 101*GeV ) { +++mumuZCR; + _hist_St_mumu_ZCR->fill(St_ll, weight); + } + } + // eejj, Z control region + else if ( cand_e.size() == 2 ) { + if ( M_ll >= 81*GeV && M_ll <= 101*GeV ) { +++eeZCR; + _hist_St_ee_ZCR->fill(St_ll, weight); + + } + } + + if ( cand_mu.size() == 1 ) { + // munujj, W+2jets control region + if ( recon_jets.size() == 2 && + mT >= 40*GeV && mT <= 150*GeV ) { +++munuW2CR; + _hist_MLQ_munu_W2CR->fill(M_LQ, weight); + } + // munujj, tt control region + if ( recon_jets.size() >= 4 && + recon_jets[0].momentum().pT() > 50*GeV && recon_jets[1].momentum().pT() > 40*GeV && recon_jets[2].momentum().pT() > 30*GeV ) { +++munuttCR; + _hist_MLQ_munu_ttCR->fill(M_LQ, weight); + } + } + if ( cand_e.size() == 1 ) { + // enujj, W+2jets control region + if ( recon_jets.size() == 2 && + mT >= 40*GeV && mT <= 150*GeV ) { +++enuW2CR; + _hist_MLQ_enu_W2CR->fill(M_LQ, weight); + } + // enujj, tt control region + if ( recon_jets.size() >= 4 && + recon_jets[0].momentum().pT() > 50*GeV && recon_jets[1].momentum().pT() > 40*GeV && recon_jets[2].momentum().pT() > 30*GeV ) { +++enuttCR; + _hist_MLQ_enu_ttCR->fill(M_LQ, weight); + } + } + + + + + // ========= PRESELECTION ======================= + + + + // Single lepton channel cuts + if ( single_lept ) { + + if ( eTmiss <= 25*GeV ) { +cerr << " ->Event vetoed. eTmiss=" << eTmiss << '\n'; + vetoEvent; + } +++eTmisscut; + + if ( mT <= 40*GeV ) + vetoEvent; + +//++mTcut; + + // enujj channel + if ( cand_e.size() == 1 && cand_mu.empty() ) { + + // Triangle cut + double dPhi_jet1 = deltaPhi( recon_jets[0].phi(), pTmiss.phi() ); + double dPhi_jet2 = deltaPhi( recon_jets[1].phi(), pTmiss.phi() ); + + if ( dPhi_jet1 <= 1.5 * (1 - eTmiss/45) || + dPhi_jet2 <= 1.5 * (1 - eTmiss/45) ) { +++emuvjj; + vetoEvent; + } + } + } + + // ==================== FILL ==================== + + + // mumujj channel + if ( cand_mu.size() == 2 ) { + if ( M_ll <= 120*GeV || + M_LQ <= 150*GeV || + p_l1.pT() <= 30*GeV || p_l2.pT() <= 30*GeV || + p_j[0].pT() <= 30*GeV || p_j[1].pT() <= 30*GeV || + St_ll <= 450*GeV ) { +cerr<<" ->Dilept event vetoed. Table 4 cuts." << '\n'; + vetoEvent; + } + else { + + +++mumujj; +cerr<< " ->MUMUJJ event selected." << '\n'; + _hist_St_mumu->fill(St_ll, weight); + _count_mumujj->fill(0.5, weight); + + } + } + // eejj channel + else if ( cand_e.size() == 2 ) { + if ( M_ll <= 120*GeV || + M_LQ <= 150*GeV || + p_l1.pT() <= 30*GeV || p_l2.pT() <= 30*GeV || + p_j[0].pT() <= 30*GeV || p_j[1].pT() <= 30*GeV || + St_ll <= 450*GeV ) { +cerr<<" ->Dilept event vetoed. Table 4 cuts." << '\n'; + vetoEvent; + } + else { + +++eejj; +cerr<< " ->EEJJ event selected." << '\n'; + _hist_St_ee->fill(St_ll, weight); + _count_eejj->fill(0.5, weight); + + } + } + + + // muvjj channel + else if ( cand_mu.size() == 1 ) { + + + + if (M_LQ<=150*GeV) { +cerr<<" ->muvjj event vetoed. Not enough M_LQ: " << M_LQ<< '\n'; + vetoEvent; + } +++MLQonelept; + if (Mt_LQ<=150*GeV) { +cerr<<" ->muvjj event vetoed. Not enough Mt_LQ: " << Mt_LQ<< '\n'; + vetoEvent; + } +++MtLQonelept; + if (St_v<=400*GeV) { +cerr<<" ->muvjj event vetoed. Not enough St_v: " << St_v<< '\n'; + vetoEvent; + } +++Stvonelept; + if (mT<=160*GeV) { +cerr<<" ->muvjj event vetoed. Not enough mT: " << mT<<'\n'; + vetoEvent; + } +++mTonelept; + //else { +++muvjj; +cerr<< " ->MUVJJ event selected." << '\n'; + _hist_MLQ_muv->fill(M_LQ, weight); + _count_muvjj->fill(0.5, weight); + + //} + } + + // evjj channel + else if ( cand_e.size() == 1 ) { + +if (M_LQ<=180*GeV) { +cerr<<" ->evjj event vetoed. Not enough M_LQ: " << M_LQ<< '\n'; + vetoEvent; + } +++MLQev; + if (Mt_LQ<=180*GeV) { +cerr<<" ->evjj event vetoed. Not enough Mt_LQ: " << Mt_LQ<< '\n'; + vetoEvent; + } +++MtLQev; + if (St_v<=410*GeV) { +cerr<<" ->evjj event vetoed. Not enough St_v: " << St_v<< '\n'; + vetoEvent; + } +++Stvev; +if (mT<=200*GeV) { +cerr<<" ->evjj event vetoed. Not enough mT: " << mT<<'\n'; + vetoEvent; + } +++mTev; + //else { +++evjj; +cerr<< " ->EVJJ event selected." << '\n'; +_hist_MLQ_ev->fill(M_LQ, weight); + _count_evjj->fill(0.5, weight); + + +/* + if ( mT <= 200*GeV || + M_LQ <= 180*GeV || + Mt_LQ <= 180*GeV || + St_v <= 410*GeV ) { +cerr<<" ->evjj event vetoed. Doesn't pass table 4 cuts." << '\n'; + vetoEvent; + } + else { +++evjj; +cerr<< " ->EVJJ event selected." << '\n'; +_hist_MLQ_ev->fill(M_LQ, weight); + _count_evjj->fill(0.5, weight); + + } + +*/ + } + + + } + + //@} + + + void finalize() { +cerr << '\n' << "Of " << count << " events, saw " +<< vetoe << " (after veto region cut), " +<< Njetscut << " (after 2jet req). " +<< '\n' +<< "For " << dilept << " dilept events: " +<< candmumujj << " cand mumujj events, " +<< candeejj << " cand eejj events." +<< '\n' +<< "For " << onelept << " onelept events: " +<< candmvjj << " preselected mvjj events, " +<< candevjj << " preselected evjj events; " +<< eTmisscut << " (eTmiss req); " +<< emuvjj << " leftover; " +<< MLQonelept << " (muvjj M_LQ cut), " +<< MtLQonelept << " (muvjj Mt_LQ cut), " +<< Stvonelept << " (muvjj St_v cut), " +<< mTonelept << " (muvjj mT cut); " +<< MLQev << " (evjj M_LQ cut), " +<< MtLQev << " (evjj Mt_LQ cut), " +<< Stvev << " (evjj St_v cut), " +<< mTev << " (evjj mT cut). " +<< '\n'<<'\n' +; + +cerr << "CR - " << "mumu Z: " << mumuZCR << " ee Z: " << eeZCR << " munu W+2jets: " << munuW2CR << " munu tt: " << munuttCR << " enu W+2jets: " << enuW2CR << " enu tt: " << enuttCR << '\n'; + +cerr << "mumujj: " << mumujj << " eejj: " << eejj << " muvjj: " << muvjj << " evjj: " << evjj << '\n'; + + + scale( _hist_St_ee, 120. * 35. * crossSection()/sumOfWeights() ); + scale( _hist_St_mumu, 120. * 35. * crossSection()/sumOfWeights() ); + scale( _hist_MLQ_muv, 50. * 35. * crossSection()/sumOfWeights() ); + scale( _hist_MLQ_ev, 50. * 35. * crossSection()/sumOfWeights() ); + + + + scale( _hist_St_mumu_ZCR, 20. * 35. * crossSection()/sumOfWeights() ); + scale( _hist_St_ee_ZCR, 20. * 35. * crossSection()/sumOfWeights() ); + scale( _hist_MLQ_munu_W2CR, 20. * 35. * crossSection()/sumOfWeights() ); + scale( _hist_MLQ_enu_W2CR, 20. * 35. * crossSection()/sumOfWeights() ); + scale( _hist_MLQ_munu_ttCR, 20. * 35. * crossSection()/sumOfWeights() ); + scale( _hist_MLQ_enu_ttCR, 20. * 35. * crossSection()/sumOfWeights() ); + +/* +scale( _hist_eTmiss_mu, binwidth*luminosity* crossSection()/sumOfWeights() ); +*/ + + } + + private: + + /// @name Histograms + //@{ + AIDA::IHistogram1D* _count_mumujj; + AIDA::IHistogram1D* _count_eejj; + AIDA::IHistogram1D* _count_muvjj; + AIDA::IHistogram1D* _count_evjj; + + AIDA::IHistogram1D* _hist_St_mumu; + AIDA::IHistogram1D* _hist_St_ee; + AIDA::IHistogram1D* _hist_MLQ_muv; + AIDA::IHistogram1D* _hist_MLQ_ev; + + AIDA::IHistogram1D* _hist_St_mumu_ZCR; + AIDA::IHistogram1D* _hist_St_ee_ZCR; + AIDA::IHistogram1D* _hist_MLQ_munu_W2CR; + AIDA::IHistogram1D* _hist_MLQ_enu_W2CR; + AIDA::IHistogram1D* _hist_MLQ_munu_ttCR; + AIDA::IHistogram1D* _hist_MLQ_enu_ttCR; + + + + + //@} + + +//DEBUG VARIABLES + +int count; +int vetoe; +int Njetscut; +int dilept; +int candmumujj; +int candeejj; +int onelept; +int eTmisscut; +int candmvjj; +int candevjj; +int mumujj; +int eejj; +int mTonelept; +int MLQonelept; +int MtLQonelept; +int Stvonelept; +int mTev; +int MLQev; +int MtLQev; +int Stvev; +int muvjj; +int evjj; +int emuvjj; +int cande; +int candmu; +int tmpe; +int tmpmu; +int mumuZCR; +int eeZCR; +int munuW2CR; +int munuttCR; +int enuW2CR; +int enuttCR; + + }; + + + + // This global object acts as a hook for the plugin system + AnalysisBuilder<ATLAS_2011_S9041966> plugin_ATLAS_2011_S9041966; + + +} Modified: trunk/src/Analyses/Makefile.am ============================================================================== --- trunk/src/Analyses/Makefile.am Tue Aug 2 13:53:41 2011 (r3274) +++ trunk/src/Analyses/Makefile.am Tue Aug 2 16:27:14 2011 (r3275) @@ -59,7 +59,9 @@ RivetATLASAnalyses_la_SOURCES += \ ATLAS_2010_CONF_2010_049.cc \ ATLAS_2011_S9019561.cc \ - ATLAS_2011_CONF_2011_090.cc + ATLAS_2011_S9041966.cc \ + ATLAS_2011_CONF_2011_090.cc \ + ATLAS_2011_CONF_2011_098.cc endif
More information about the Rivet-svn mailing list |