rivet is hosted by Hepforge, IPPP Durham
ATLAS_2012_I1186556.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 #include "Rivet/Analysis.hh"
00003 #include "Rivet/Tools/BinnedHistogram.hh"
00004 #include "Rivet/RivetYODA.hh"
00005 #include "Rivet/Tools/Logging.hh"
00006 #include "Rivet/Projections/FinalState.hh"
00007 #include "Rivet/Projections/ChargedFinalState.hh"
00008 #include "Rivet/Projections/VisibleFinalState.hh"
00009 #include "Rivet/Projections/IdentifiedFinalState.hh"
00010 #include "Rivet/Projections/VetoedFinalState.hh"
00011 #include "Rivet/Projections/FastJets.hh"
00012 #include "Rivet/Tools/RivetMT2.hh"
00013 
00014 namespace Rivet {
00015 
00016 
00017   class ATLAS_2012_I1186556 : public Analysis {
00018   public:
00019 
00020     /// @name Constructors etc.
00021     //@{
00022 
00023     /// Constructor
00024 
00025     ATLAS_2012_I1186556()
00026       : Analysis("ATLAS_2012_I1186556")
00027     {    }
00028 
00029     //@}
00030 
00031 
00032   public:
00033 
00034     /// @name Analysis methods
00035     //@{
00036 
00037     /// Book histograms and initialize projections before the run
00038     void init() {
00039 
00040       // projection to find the electrons
00041       vector<pair<double, double> > eta_e;
00042       eta_e.push_back(make_pair(-2.47,2.47));
00043       IdentifiedFinalState elecs(eta_e, 20.0*GeV);
00044       elecs.acceptIdPair(ELECTRON);
00045       addProjection(elecs, "elecs");
00046 
00047       // projection to find the muons
00048       vector<pair<double, double> > eta_m;
00049       eta_m.push_back(make_pair(-2.4,2.4));
00050       IdentifiedFinalState muons(eta_m, 10.0*GeV);
00051       muons.acceptIdPair(MUON);
00052       addProjection(muons, "muons");
00053 
00054       // Jet finder
00055       VetoedFinalState vfs;
00056       vfs.addVetoPairId(MUON);
00057       addProjection(FastJets(vfs, FastJets::ANTIKT, 0.4), "AntiKtJets04");
00058 
00059       // all tracks (to do deltaR with leptons)
00060       addProjection(ChargedFinalState(-3.0,3.0,1.*GeV),"cfs");
00061 
00062       // for pTmiss
00063       addProjection(VisibleFinalState(-4.9,4.9),"vfs");
00064 
00065       // Book histograms
00066       _count_SR_SF     = bookHisto1D("count_SR_SF"    , 1, 0., 1.);
00067       _count_SR_OF     = bookHisto1D("count_SR_OF"    , 1, 0., 1.);
00068 
00069       _hist_mT2_SF_exp = bookHisto1D("hist_mT2_SF_exp", 40 , 0., 200. );
00070       _hist_mT2_OF_exp = bookHisto1D("hist_mT2_OF_exp", 40 , 0., 200. );
00071       _hist_mT2_SF_MC  = bookHisto1D("hist_mT2_SF_MC" , 500, 0., 1000.);
00072       _hist_mT2_OF_MC  = bookHisto1D("hist_mT2_OF_MC" , 500, 0., 1000.);
00073 
00074     }
00075 
00076     /// Perform the per-event analysis
00077     void analyze(const Event& event) {
00078       const double weight = event.weight();
00079 
00080       // get the candiate jets
00081       Jets cand_jets;
00082       foreach ( const Jet& jet,
00083                 applyProjection<FastJets>(event, "AntiKtJets04").jetsByPt(20.0*GeV) ) {
00084         if ( fabs( jet.momentum().eta() ) < 4.5 ) {
00085           cand_jets.push_back(jet);
00086         }
00087       }
00088       // charged tracks for isolation
00089       ParticleVector chg_tracks =
00090         applyProjection<ChargedFinalState>(event, "cfs").particles();
00091       // find the electrons
00092       ParticleVector cand_e;
00093       foreach( const Particle & e,
00094                applyProjection<IdentifiedFinalState>(event, "elecs").particlesByPt()) {
00095         // remove any leptons within 0.4 of any candidate jets
00096         bool e_near_jet = false;
00097         foreach ( const Jet& jet, cand_jets ) {
00098           double dR = deltaR(e.momentum(),jet.momentum());
00099           if ( dR < 0.4 && dR > 0.2 ) {
00100             e_near_jet = true;
00101             break;
00102           }
00103         }
00104     if ( e_near_jet ) continue;
00105     cand_e.push_back(e);
00106       }
00107       ParticleVector cand_mu;
00108       foreach( const Particle & mu,
00109                applyProjection<IdentifiedFinalState>(event, "muons").particlesByPt()) {
00110         // remove any leptons within 0.4 of any candidate jets
00111         bool mu_near_jet = false;
00112         foreach ( const Jet& jet, cand_jets ) {
00113           if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 ) {
00114             mu_near_jet = true;
00115             break;
00116           }
00117         }
00118         if ( mu_near_jet ) continue;
00119     cand_mu.push_back(mu);
00120       }
00121       // pTcone around muon track
00122       ParticleVector recon_mu;
00123       foreach ( const Particle & mu, cand_mu ) {
00124         double pTinCone = -mu.momentum().pT();
00125         foreach ( const Particle & track, chg_tracks ) {
00126           if ( deltaR(mu.momentum(),track.momentum()) < 0.2 )
00127             pTinCone += track.momentum().pT();
00128         }
00129         if ( pTinCone < 1.8*GeV ) recon_mu.push_back(mu);
00130       }
00131       // pTcone around electron track
00132       ParticleVector recon_e;
00133       foreach ( const Particle & e, cand_e ) {
00134         double pTinCone = -e.momentum().pT();
00135         foreach ( const Particle & track, chg_tracks ) {
00136           if ( deltaR(e.momentum(),track.momentum()) < 0.2 )
00137             pTinCone += track.momentum().pT();
00138         }
00139         if ( pTinCone < 0.1 * e.momentum().pT() ) recon_e.push_back(e);
00140       }
00141 
00142       // pTmiss
00143       FourMomentum pTmiss;
00144       foreach ( const Particle & p,
00145                 applyProjection<VisibleFinalState>(event, "vfs").particles() ) {
00146         pTmiss -= p.momentum();
00147       }
00148 
00149       // discard jets that overlap with electrons
00150       Jets recon_jets;
00151       foreach ( const Jet& jet, cand_jets ) {
00152         if(fabs(jet.momentum().eta())>2.5||
00153            jet.momentum().perp()<20.) continue;
00154     bool away_from_e = true;
00155     foreach ( const Particle & e, cand_e ) {
00156       if ( deltaR(e.momentum(),jet.momentum()) < 0.2 ) {
00157         away_from_e = false;
00158         break;
00159       }
00160     }
00161     if ( away_from_e ) recon_jets.push_back( jet );
00162       }
00163 
00164       // put leptons into 1 vector and order by pT
00165       ParticleVector leptons(recon_e.begin(),recon_e.end());
00166       leptons.insert(leptons.begin(),recon_mu.begin(),recon_mu.end());
00167       sort(leptons.begin(),leptons.end(),cmpParticleByPt);
00168 
00169       // exactly two leptons
00170       if(leptons.size() !=2) vetoEvent;
00171 
00172       // hardest lepton pT greater the 25 (20) e(mu)
00173       if( (abs(leptons[0].pdgId())==ELECTRON && leptons[0].momentum().perp()<25.) ||
00174       (abs(leptons[0].pdgId())==ELECTRON && leptons[0].momentum().perp()<20.))
00175     vetoEvent;
00176 
00177       // require opposite sign
00178       if(leptons[0].pdgId()*leptons[1].pdgId()>0) vetoEvent;
00179 
00180       // and invariant mass > 20
00181       double mll = (leptons[0].momentum()+leptons[1].momentum()).mass();
00182       if(mll<20.) vetoEvent;
00183 
00184       // two jets 1st pT > 50 and second pT> 25
00185       if(recon_jets.size()<2 || recon_jets[0].momentum().perp()<50. ||
00186      recon_jets[1].momentum().perp()<25.) vetoEvent;
00187 
00188       // calculate mT2
00189       double m_T2 = mT2::mT2( leptons[0].momentum(),leptons[1].momentum(),
00190                   pTmiss,0.0 ); // zero mass invisibles
00191 
00192       // same flavour region
00193       if(leptons[0].pdgId()==-leptons[1].pdgId()) {
00194     // remove Z region
00195     if(mll>71.&&mll<111.) vetoEvent;
00196     // require at least 1 b jet
00197     unsigned int n_b=0;
00198     for(unsigned int ix=0;ix<recon_jets.size();++ix) {
00199        if(recon_jets[ix].containsBottom() && rand()/static_cast<double>(RAND_MAX)<=0.60)
00200          ++n_b;
00201     }
00202     if(n_b==0) vetoEvent;
00203     _hist_mT2_SF_exp->fill(m_T2,weight);
00204     _hist_mT2_SF_MC ->fill(m_T2,weight);
00205     if(m_T2>120.) _count_SR_SF->fill(0.5,weight);
00206       }
00207       // opposite flavour region
00208       else {
00209     _hist_mT2_OF_exp->fill(m_T2,weight);
00210     _hist_mT2_OF_MC ->fill(m_T2,weight);
00211     if(m_T2>120.) _count_SR_OF->fill(0.5,weight);
00212       }
00213     }
00214     //@}
00215 
00216 
00217     void finalize() {
00218 
00219       double norm = 4.7* crossSection()/sumOfWeights()/femtobarn;
00220       scale(_count_SR_SF    ,   norm);
00221       scale(_count_SR_OF    ,   norm);
00222       scale(_hist_mT2_SF_exp,5.*norm);
00223       scale(_hist_mT2_OF_exp,5.*norm);
00224       scale(_hist_mT2_SF_MC ,   norm/4.7);
00225       scale(_hist_mT2_OF_MC ,   norm/4.7);
00226 
00227     }
00228 
00229   private:
00230 
00231     /// @name Histograms
00232     //@{
00233     Histo1DPtr _count_SR_SF;
00234     Histo1DPtr _count_SR_OF;
00235 
00236     Histo1DPtr _hist_mT2_SF_exp;
00237     Histo1DPtr _hist_mT2_OF_exp;
00238     Histo1DPtr _hist_mT2_SF_MC;
00239     Histo1DPtr _hist_mT2_OF_MC;
00240     //@}
00241 
00242   };
00243 
00244   // The hook for the plugin system
00245   DECLARE_RIVET_PLUGIN(ATLAS_2012_I1186556);
00246 
00247 }