diff --git a/analyses/pluginATLAS/ATLAS_2014_I1298023.cc b/analyses/pluginATLAS/ATLAS_2014_I1298023.cc --- a/analyses/pluginATLAS/ATLAS_2014_I1298023.cc +++ b/analyses/pluginATLAS/ATLAS_2014_I1298023.cc @@ -1,128 +1,126 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/ChargedLeptons.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/MissingMomentum.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/DressedLeptons.hh" namespace Rivet { /// ATLAS same-sign WW production at 8 TeV (PRL) class ATLAS_2014_I1298023 : public Analysis { public: /// @name Constructors etc. //@{ ATLAS_2014_I1298023() : Analysis("ATLAS_2014_I1298023") { setNeedsCrossSection(true); } //@} public: /// @name Analysis methods //@{ void init() { const FinalState fs; // bare leptons ChargedLeptons bare_leptons(fs); // dressed leptons Cut cuts = (Cuts::abseta < 2.5) & (Cuts::pT > 25*GeV); DressedLeptons leptons(fs, bare_leptons, 0.1, cuts); declare(leptons, "leptons"); // MET declare(MissingMomentum(fs), "MissingET"); // jets VetoedFinalState vfs(fs); vfs.addVetoPairId(PID::MUON); vfs.vetoNeutrinos(); declare(FastJets(vfs, FastJets::ANTIKT, 0.4), "jets"); // book histogram book(_hist ,1, 1, 1); } /// Perform the per-event analysis void analyze(const Event& event) { - const double weight = 1.0; - const vector& leptons = apply(event, "leptons").dressedLeptons(); if ( leptons.size() < 2 ) vetoEvent; double minDR_ll = MAXDOUBLE, mll = -1.0; for (unsigned int i = 0; i < leptons.size(); ++i) { DressedLepton lep1 = leptons[i]; for (unsigned int j = i + 1; j < leptons.size(); ++j) { DressedLepton lep2 = leptons[j]; double dr = deltaR(lep1, lep2); if ( dr < minDR_ll ) minDR_ll = dr; double m = (lep1.momentum() + lep2.momentum()).mass(); if ( mll < 0. || m < mll ) mll = m; } } if ( minDR_ll <= 0.3 || mll <= 20*GeV ) vetoEvent; if ( leptons[0].charge() * leptons[1].charge() < 0.0 ) vetoEvent; const MissingMomentum& met = apply(event, "MissingET"); if ( met.vectorEt().mod() <= 40*GeV ) vetoEvent; const Jets& all_jets = apply(event, "jets").jetsByPt( (Cuts::abseta < 4.5) && (Cuts::pT > 30*GeV) ); Jets jets; double minDR_overall = MAXDOUBLE; foreach (const Jet& jet, all_jets) { double minDR_jet = MAXDOUBLE, minDR_electrons = MAXDOUBLE; foreach( DressedLepton lep, leptons ) { double dr = deltaR(jet, lep); if ( dr < minDR_jet ) minDR_jet = dr; if ( lep.abspid() == 11 && dr < minDR_electrons ) minDR_electrons = dr; } if ( minDR_electrons < 0.05 ) continue; // veto jet if it overlaps with electron if ( minDR_jet < minDR_overall ) minDR_overall = minDR_jet; jets += jet; } if ( jets.size() < 2 || minDR_overall <= 0.3 ) vetoEvent; FourMomentum dijet = jets[0].momentum() + jets[1].momentum(); if ( dijet.mass() <= 500*GeV ) vetoEvent; // inclusive region - _hist->fill(0.5, weight); + _hist->fill(0.5); // VBS region - if ( deltaRap(jets[0], jets[1]) > 2.4 ) _hist->fill(1.5, weight); + if ( deltaRap(jets[0], jets[1]) > 2.4 ) _hist->fill(1.5); } /// Normalise histograms etc., after the run void finalize() { const double scalefactor( crossSection() / sumOfWeights() ); scale(_hist, scalefactor); } //@} private: /// @name Histograms //@{ Histo1DPtr _hist; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2014_I1298023); } diff --git a/analyses/pluginATLAS/ATLAS_2014_I1298811.cc b/analyses/pluginATLAS/ATLAS_2014_I1298811.cc --- a/analyses/pluginATLAS/ATLAS_2014_I1298811.cc +++ b/analyses/pluginATLAS/ATLAS_2014_I1298811.cc @@ -1,199 +1,196 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { class ATLAS_2014_I1298811 : public Analysis { public: ATLAS_2014_I1298811() : Analysis("ATLAS_2014_I1298811") { } void init() { // Configure projections const FinalState fs(-4.8, 4.8, 0*MeV); declare(fs, "FS"); const FastJets jets(fs, FastJets::ANTIKT, 0.4); declare(jets, "Jets"); // Book histograms for (size_t itopo = 0; itopo < 2; ++itopo) { // Profiles for (size_t iregion = 0; iregion < 3; ++iregion) { book(_p_ptsumch_vs_ptlead[itopo][iregion] ,1+iregion, 1, itopo+1); book(_p_nch_vs_ptlead[itopo][iregion] ,4+iregion, 1, itopo+1); } book(_p_etsum25_vs_ptlead_trans[itopo] ,7, 1, itopo+1); book(_p_etsum48_vs_ptlead_trans[itopo] ,8, 1, itopo+1); book(_p_chratio_vs_ptlead_trans[itopo] ,9, 1, itopo+1); book(_p_ptmeanch_vs_ptlead_trans[itopo] ,10, 1, itopo+1); book(_p_ptmeanch_vs_nch_trans[0] ,11, 1, 1); book(_p_ptmeanch_vs_nch_trans[1] ,12, 1, 1); // 1D histos for (size_t iregion = 0; iregion < 3; ++iregion) { for (size_t ipt = 0; ipt < 4; ++ipt) { book(_h_ptsumch[ipt][itopo][iregion] ,13+3*ipt+iregion, 1, itopo+1); book(_h_nch[ipt][itopo][iregion] ,25+3*ipt+iregion, 1, itopo+1); } } } } void analyze(const Event& event) { // Find the jets with pT > 20 GeV and *rapidity* within 2.8 /// @todo Use Cuts instead rather than an eta cut in the proj and a y cut after const Jets alljets = apply(event, "Jets").jetsByPt(20*GeV); Jets jets; foreach (const Jet& j, alljets) if (j.absrap() < 2.8) jets.push_back(j); // Require at least one jet in the event if (jets.empty()) vetoEvent; - // Get the event weight since we will be filling some histos - const double weight = 1.0; - // Identify the leading jet and its phi and pT const FourMomentum plead = jets[0].momentum(); const double philead = plead.phi(); const double etalead = plead.eta(); const double ptlead = plead.pT(); MSG_DEBUG("Leading object: pT = " << ptlead << ", eta = " << etalead << ", phi = " << philead); // Sum particle properties in the transverse regions int tmpnch[2] = {0,0}; double tmpptsum[2] = {0,0}; double tmpetsum48[2] = {0,0}; double tmpetsum25[2] = {0,0}; const Particles particles = apply(event, "FS").particles(); foreach (const Particle& p, particles) { // Only consider the transverse region(s), not toward or away if (!inRange(deltaPhi(p.phi(), philead), PI/3.0, TWOPI/3.0)) continue; // Work out which transverse side this particle is on const size_t iside = (mapAngleMPiToPi(p.phi() - philead) > 0) ? 0 : 1; MSG_TRACE(p.phi() << " vs. " << philead << ": " << iside); // Charged or neutral particle? const bool charged = PID::threeCharge(p.pdgId()) != 0; // Track observables if (charged && fabs(p.eta()) < 2.5 && p.pT() > 500*MeV) { tmpnch[iside] += 1; tmpptsum[iside] += p.pT(); } // Cluster observables if ((charged && p.p3().mod() > 200*MeV) || (!charged && p.p3().mod() > 500*MeV)) { tmpetsum48[iside] += p.pT(); if (fabs(p.eta()) < 2.5) tmpetsum25[iside] += p.pT(); } } // Construct tot/max/min counts (for trans/max/min, indexed by iregion) const int nch[3] = { tmpnch[0] + tmpnch[1], std::max(tmpnch[0], tmpnch[1]), std::min(tmpnch[0], tmpnch[1]) }; const double ptsum[3] = { tmpptsum[0] + tmpptsum[1], std::max(tmpptsum[0], tmpptsum[1]), std::min(tmpptsum[0], tmpptsum[1]) }; const double etsum48[3] = { tmpetsum48[0] + tmpetsum48[1], std::max(tmpetsum48[0], tmpetsum48[1]), std::min(tmpetsum48[0], tmpetsum48[1]) }; const double etsum25[3] = { tmpetsum25[0] + tmpetsum25[1], std::max(tmpetsum25[0], tmpetsum25[1]), std::min(tmpetsum25[0], tmpetsum25[1]) }; ////////////////////////////////////////////////////////// // Now fill the histograms with the computed quantities // phi sizes of each trans/max/min region (for indexing by iregion) const double dphi[3] = { 2*PI/3.0, PI/3.0, PI/3.0 }; // Loop over inclusive jet and exclusive dijet configurations for (size_t itopo = 0; itopo < 2; ++itopo) { // Exit early if in the exclusive dijet iteration and the exclusive dijet cuts are not met if (itopo == 1) { if (jets.size() != 2) continue; const FourMomentum psublead = jets[1].momentum(); // Delta(phi) cut const double phisublead = psublead.phi(); if (deltaPhi(philead, phisublead) < 2.5) continue; // pT fraction cut const double ptsublead = psublead.pT(); if (ptsublead < 0.5*ptlead) continue; MSG_DEBUG("Exclusive dijet event"); } // Plot profiles and distributions which have no max/min region definition - _p_etsum25_vs_ptlead_trans[itopo]->fill(ptlead/GeV, etsum25[0]/5.0/dphi[0]/GeV, weight); - _p_etsum48_vs_ptlead_trans[itopo]->fill(ptlead/GeV, etsum48[0]/9.6/dphi[0]/GeV, weight); + _p_etsum25_vs_ptlead_trans[itopo]->fill(ptlead/GeV, etsum25[0]/5.0/dphi[0]/GeV); + _p_etsum48_vs_ptlead_trans[itopo]->fill(ptlead/GeV, etsum48[0]/9.6/dphi[0]/GeV); if (etsum25[0] > 0) { - _p_chratio_vs_ptlead_trans[itopo]->fill(ptlead/GeV, ptsum[0]/etsum25[0], weight); + _p_chratio_vs_ptlead_trans[itopo]->fill(ptlead/GeV, ptsum[0]/etsum25[0]); } const double ptmean = safediv(ptsum[0], nch[0], -1); ///< Return -1 if div by zero if (ptmean >= 0) { - _p_ptmeanch_vs_ptlead_trans[itopo]->fill(ptlead/GeV, ptmean/GeV, weight); - _p_ptmeanch_vs_nch_trans[itopo]->fill(nch[0], ptmean/GeV, weight); + _p_ptmeanch_vs_ptlead_trans[itopo]->fill(ptlead/GeV, ptmean/GeV); + _p_ptmeanch_vs_nch_trans[itopo]->fill(nch[0], ptmean/GeV); } // Plot remaining profile and 1D observables, which are defined in all 3 tot/max/min regions for (size_t iregion = 0; iregion < 3; ++iregion) { - _p_ptsumch_vs_ptlead[itopo][iregion]->fill(ptlead/GeV, ptsum[iregion]/5.0/dphi[iregion]/GeV, weight); - _p_nch_vs_ptlead[itopo][iregion]->fill(ptlead/GeV, nch[iregion]/5.0/dphi[iregion], weight); + _p_ptsumch_vs_ptlead[itopo][iregion]->fill(ptlead/GeV, ptsum[iregion]/5.0/dphi[iregion]/GeV); + _p_nch_vs_ptlead[itopo][iregion]->fill(ptlead/GeV, nch[iregion]/5.0/dphi[iregion]); for (size_t ipt = 0; ipt < 4; ++ipt) { if (ipt == 1 && !inRange(ptlead/GeV, 20, 60)) continue; if (ipt == 2 && !inRange(ptlead/GeV, 60, 210)) continue; if (ipt == 3 && ptlead/GeV < 210) continue; - _h_ptsumch[ipt][itopo][iregion]->fill(ptsum[iregion]/5.0/dphi[iregion]/GeV, weight); - _h_nch[ipt][itopo][iregion]->fill(nch[iregion]/5.0/dphi[iregion], weight); + _h_ptsumch[ipt][itopo][iregion]->fill(ptsum[iregion]/5.0/dphi[iregion]/GeV); + _h_nch[ipt][itopo][iregion]->fill(nch[iregion]/5.0/dphi[iregion]); } } } } void finalize() { for (size_t iregion = 0; iregion < 3; ++iregion) { for (size_t itopo = 0; itopo < 2; ++itopo) { for (size_t ipt = 0; ipt < 4; ++ipt) { normalize(_h_ptsumch[ipt][itopo][iregion], 1.0); normalize(_h_nch[ipt][itopo][iregion], 1.0); } } } } private: /// @name Histogram arrays //@{ Profile1DPtr _p_ptsumch_vs_ptlead[2][3]; Profile1DPtr _p_nch_vs_ptlead[2][3]; Profile1DPtr _p_ptmeanch_vs_ptlead_trans[2]; Profile1DPtr _p_etsum25_vs_ptlead_trans[2]; Profile1DPtr _p_etsum48_vs_ptlead_trans[2]; Profile1DPtr _p_chratio_vs_ptlead_trans[2]; Profile1DPtr _p_ptmeanch_vs_nch_trans[2]; Histo1DPtr _h_ptsumch[4][2][3]; Histo1DPtr _h_nch[4][2][3]; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2014_I1298811); } diff --git a/analyses/pluginATLAS/ATLAS_2014_I1304688.cc b/analyses/pluginATLAS/ATLAS_2014_I1304688.cc --- a/analyses/pluginATLAS/ATLAS_2014_I1304688.cc +++ b/analyses/pluginATLAS/ATLAS_2014_I1304688.cc @@ -1,325 +1,324 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/FastJets.hh" #include namespace Rivet { /// @brief ATLAS 7 TeV jets in ttbar events analysis /// /// @author W. H. Bell /// @author A. Grohsjean class ATLAS_2014_I1304688 : public Analysis { public: ATLAS_2014_I1304688(): Analysis("ATLAS_2014_I1304688"), _jet_ntag(0), _met_et(0.), _met_phi(0.), _hMap(), //_chanLimit(3), _histLimit(6) { } void init() { // Eta ranges /// @todo 1 MeV? Really? Cut eta_full = Cuts::abseta < 5.0 && Cuts::pT > 1.0*MeV; Cut eta_lep = Cuts::abseta < 2.5; // All final state particles FinalState fs(eta_full); // Get photons to dress leptons IdentifiedFinalState photons(fs); photons.acceptIdPair(PID::PHOTON); // Projection to find the electrons IdentifiedFinalState el_id(fs); el_id.acceptIdPair(PID::ELECTRON); PromptFinalState electrons(el_id); electrons.acceptTauDecays(true); declare(electrons, "electrons"); DressedLeptons dressedelectrons(photons, electrons, 0.1, eta_lep && Cuts::pT > 25*GeV, true); declare(dressedelectrons, "dressedelectrons"); DressedLeptons vetodressedelectrons(photons, electrons, 0.1, eta_lep && Cuts::pT >= 15*GeV, true); declare(vetodressedelectrons, "vetodressedelectrons"); DressedLeptons ewdressedelectrons(photons, electrons, 0.1, eta_full, true); declare(ewdressedelectrons, "ewdressedelectrons"); // Projection to find the muons IdentifiedFinalState mu_id(fs); mu_id.acceptIdPair(PID::MUON); PromptFinalState muons(mu_id); muons.acceptTauDecays(true); declare(muons, "muons"); vector > eta_muon; DressedLeptons dressedmuons(photons, muons, 0.1, eta_lep && Cuts::pT >= 25*GeV, true); declare(dressedmuons, "dressedmuons"); DressedLeptons vetodressedmuons(photons, muons, 0.1, eta_lep && Cuts::pT >= 15*GeV, true); declare(vetodressedmuons, "vetodressedmuons"); DressedLeptons ewdressedmuons(photons, muons, 0.1, eta_full, true); declare(ewdressedmuons, "ewdressedmuons"); // Projection to find neutrinos and produce MET IdentifiedFinalState nu_id; nu_id.acceptNeutrinos(); PromptFinalState neutrinos(nu_id); neutrinos.acceptTauDecays(true); declare(neutrinos, "neutrinos"); // Jet clustering. VetoedFinalState vfs; vfs.addVetoOnThisFinalState(ewdressedelectrons); vfs.addVetoOnThisFinalState(ewdressedmuons); vfs.addVetoOnThisFinalState(neutrinos); FastJets jets(vfs, FastJets::ANTIKT, 0.4); jets.useInvisibles(); declare(jets, "jets"); // Book histograms for (unsigned int ihist = 0; ihist < _histLimit ; ihist++) { const unsigned int threshLimit = _thresholdLimit(ihist); for (unsigned int ithres = 0; ithres < threshLimit; ithres++) { _histogram(ihist, ithres); // Create all histograms } } } void analyze(const Event& event) { // Get the selected objects, using the projections. _dressedelectrons = sortByPt(apply(event, "dressedelectrons").dressedLeptons()); _vetodressedelectrons = apply(event, "vetodressedelectrons").dressedLeptons(); _dressedmuons = sortByPt(apply(event, "dressedmuons").dressedLeptons()); _vetodressedmuons = apply(event, "vetodressedmuons").dressedLeptons(); _neutrinos = apply(event, "neutrinos").particlesByPt(); _jets = apply(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5); // Calculate the missing ET, using the prompt neutrinos only (really?) /// @todo Why not use MissingMomentum? FourMomentum pmet; foreach (const Particle& p, _neutrinos) pmet += p.momentum(); _met_et = pmet.pT(); _met_phi = pmet.phi(); // Check overlap of jets/leptons. unsigned int i,j; _jet_ntag = 0; _overlap = false; for (i = 0; i < _jets.size(); i++) { const Jet& jet = _jets[i]; // If dR(el,jet) < 0.4 skip the event foreach (const DressedLepton& el, _dressedelectrons) { if (deltaR(jet, el) < 0.4) _overlap = true; } // If dR(mu,jet) < 0.4 skip the event foreach (const DressedLepton& mu, _dressedmuons) { if (deltaR(jet, mu) < 0.4) _overlap = true; } // If dR(jet,jet) < 0.5 skip the event for (j = 0; j < _jets.size(); j++) { const Jet& jet2 = _jets[j]; if (i == j) continue; // skip the diagonal if (deltaR(jet, jet2) < 0.5) _overlap = true; } // Count the number of b-tags if (!jet.bTags().empty()) _jet_ntag += 1; } // Evaluate basic event selection unsigned int ejets_bits = 0, mujets_bits = 0; bool pass_ejets = _ejets(ejets_bits); bool pass_mujets = _mujets(mujets_bits); // Remove events with object overlap if (_overlap) vetoEvent; // basic event selection requirements if (!pass_ejets && !pass_mujets) vetoEvent; // Check if the additional pT threshold requirements are passed bool pass_jetPt = _additionalJetCuts(); // Count the jet multiplicity for 25, 40, 60 and 80GeV unsigned int ithres, jet_n[4]; for (ithres = 0; ithres < 4; ithres++) { jet_n[ithres] = _countJets(ithres); } // Fill histograms - const double weight = 1.0; for (unsigned int ihist = 0; ihist < 6; ihist++) { if (ihist > 0 && !pass_jetPt) continue; // additional pT threshold cuts for pT plots unsigned int threshLimit = _thresholdLimit(ihist); for (ithres = 0; ithres < threshLimit; ithres++) { if (jet_n[ithres] < 3) continue; // 3 or more jets for ljets // Fill - if (ihist == 0) _histogram(ihist, ithres)->fill(jet_n[ithres], weight); // njets - else if (ihist == 1) _histogram(ihist, ithres)->fill(_jets[0].pT(), weight); // leading jet pT - else if (ihist == 2) _histogram(ihist, ithres)->fill(_jets[1].pT(), weight); // 2nd jet pT - else if (ihist == 3 && jet_n[ithres] >= 3) _histogram(ihist, ithres)->fill(_jets[2].pT(), weight); // 3rd jet pT - else if (ihist == 4 && jet_n[ithres] >= 4) _histogram(ihist, ithres)->fill(_jets[3].pT(), weight); // 4th jet pT - else if (ihist == 5 && jet_n[ithres] >= 5) _histogram(ihist, ithres)->fill(_jets[4].pT(), weight); // 5th jet pT + if (ihist == 0) _histogram(ihist, ithres)->fill(jet_n[ithres]); // njets + else if (ihist == 1) _histogram(ihist, ithres)->fill(_jets[0].pT()); // leading jet pT + else if (ihist == 2) _histogram(ihist, ithres)->fill(_jets[1].pT()); // 2nd jet pT + else if (ihist == 3 && jet_n[ithres] >= 3) _histogram(ihist, ithres)->fill(_jets[2].pT()); // 3rd jet pT + else if (ihist == 4 && jet_n[ithres] >= 4) _histogram(ihist, ithres)->fill(_jets[3].pT()); // 4th jet pT + else if (ihist == 5 && jet_n[ithres] >= 5) _histogram(ihist, ithres)->fill(_jets[4].pT()); // 5th jet pT } } } void finalize() { // Normalize to cross-section const double norm = crossSection()/sumOfWeights(); typedef map::value_type IDtoHisto1DPtr; ///< @todo Remove when C++11 allowed foreach (IDtoHisto1DPtr ihpair, _hMap) scale(ihpair.second, norm); ///< @todo Use normalize(ihpair.second, crossSection()) // Calc averages for (unsigned int ihist = 0; ihist < _histLimit ; ihist++) { unsigned int threshLimit = _thresholdLimit(ihist); for (unsigned int ithres = 0; ithres < threshLimit; ithres++) { scale(_histogram(ihist, ithres), 0.5); } } } private: /// @name Cut helper functions //@{ // Event selection functions bool _ejets(unsigned int& cutBits) { // 1. More than zero good electrons cutBits += 1; if (_dressedelectrons.size() == 0) return false; // 2. No additional electrons passing the veto selection cutBits += 1 << 1; if (_vetodressedelectrons.size() > 1) return false; // 3. No muons passing the veto selection cutBits += 1 << 2; if (_vetodressedmuons.size() > 0) return false; // 4. total neutrino pT > 30 GeV cutBits += 1 << 3; if (_met_et <= 30.0*GeV) return false; // 5. MTW > 35 GeV cutBits += 1 << 4; if (_transMass(_dressedelectrons[0].pT(), _dressedelectrons[0].phi(), _met_et, _met_phi) <= 35*GeV) return false; // 6. At least one b-tagged jet cutBits += 1 << 5; if (_jet_ntag < 1) return false; // 7. At least three good jets cutBits += 1 << 6; if (_jets.size() < 3) return false; cutBits += 1 << 7; return true; } bool _mujets(unsigned int& cutBits) { // 1. More than zero good muons cutBits += 1; if (_dressedmuons.size() == 0) return false; // 2. No additional muons passing the veto selection cutBits += 1 << 1; if (_vetodressedmuons.size() > 1) return false; // 3. No electrons passing the veto selection cutBits += 1 << 2; if (_vetodressedelectrons.size() > 0) return false; // 4. total neutrino pT > 30 GeV cutBits += 1 << 3; if (_met_et <= 30*GeV) return false; // 5. MTW > 35 GeV cutBits += 1 << 4; if (_transMass(_dressedmuons[0].pT(), _dressedmuons[0].phi(), _met_et, _met_phi) <= 35*GeV) return false; // 6. At least one b-tagged jet cutBits += 1 << 5; if (_jet_ntag < 1) return false; // 7. At least three good jets cutBits += 1 << 6; if (_jets.size() < 3) return false; cutBits += 1 << 7; return true; } bool _additionalJetCuts() { if (_jets.size() < 2) return false; if (_jets[0].pT() <= 50*GeV || _jets[1].pT() <= 35*GeV) return false; return true; } //@} /// @name Histogram helper functions //@{ unsigned int _thresholdLimit(unsigned int histId) { if (histId == 0) return 4; return 1; } Histo1DPtr _histogram(unsigned int histId, unsigned int thresholdId) { assert(histId < _histLimit); assert(thresholdId < _thresholdLimit(histId)); const unsigned int hInd = (histId == 0) ? thresholdId : (_thresholdLimit(0) + (histId-1) + thresholdId); if (_hMap.find(hInd) != _hMap.end()) return _hMap[hInd]; _hMap.insert( { hInd, Histo1DPtr() } ); if (histId == 0) book(_hMap[hInd], 1, thresholdId+1, 1); else book(_hMap[hInd], 2, histId, 1); return _hMap[hInd]; } //@} /// @name Physics object helper functions //@{ double _transMass(double ptLep, double phiLep, double met, double phiMet) { return sqrt(2.0*ptLep*met*(1 - cos(phiLep-phiMet))); } unsigned int _countJets(unsigned int ithres) { if (ithres > 4) assert(0); double pTcut[4] = {25.,40.,60.,80.}; unsigned int i, jet_n = 0; for (i = 0; i < _jets.size(); i++) { if (_jets[i].pT() > pTcut[ithres]) jet_n++; } unsigned int ncutoff[4] = {8,7,6,5}; if (jet_n > ncutoff[ithres]) jet_n = ncutoff[ithres]; return jet_n; } //@} private: /// @name Objects that are used by the event selection decisions //@{ vector _dressedelectrons; vector _vetodressedelectrons; vector _dressedmuons; vector _vetodressedmuons; Particles _neutrinos; Jets _jets; unsigned int _jet_ntag; /// @todo Why not store the whole MET FourMomentum? double _met_et, _met_phi; bool _overlap; map _hMap; //unsigned int _chanLimit; unsigned int _histLimit; //@} }; // Declare the class as a hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2014_I1304688); } diff --git a/analyses/pluginATLAS/ATLAS_2014_I1306294.cc b/analyses/pluginATLAS/ATLAS_2014_I1306294.cc --- a/analyses/pluginATLAS/ATLAS_2014_I1306294.cc +++ b/analyses/pluginATLAS/ATLAS_2014_I1306294.cc @@ -1,233 +1,230 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ZFinder.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/HeavyHadrons.hh" #include "Rivet/Projections/VetoedFinalState.hh" namespace Rivet { class ATLAS_2014_I1306294 : public Analysis { public: /// @name Constructors etc. //@{ /// Constructors - ATLAS_2014_I1306294(std::string name="ATLAS_2014_I1306294") + ATLAS_2014_I1306294(std::string name="ATLAS_2014_I1306294") : Analysis(name) { _mode = 1; setNeedsCrossSection(true); } //@} public: /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { FinalState fs; Cut cuts = Cuts::etaIn(-2.5,2.5) & (Cuts::pT > 20.0*GeV); ZFinder zfinder(fs, cuts, _mode==1? PID::ELECTRON : PID::MUON, 76.0*GeV, 106.0*GeV, 0.1, ZFinder::CLUSTERNODECAY, ZFinder::NOTRACK); declare(zfinder, "ZFinder"); //FastJets jetpro1( getProjection("ZFinder").remainingFinalState(), FastJets::ANTIKT, 0.4); VetoedFinalState jet_fs(fs); jet_fs.addVetoOnThisFinalState(getProjection("ZFinder")); FastJets jetpro1(jet_fs, FastJets::ANTIKT, 0.4); jetpro1.useInvisibles(); declare(jetpro1, "AntiKtJets04"); declare(HeavyHadrons(), "BHadrons"); //Histograms with data binning book(_h_bjet_Pt , 3, 1, 1); book(_h_bjet_Y , 5, 1, 1); book(_h_bjet_Yboost , 7, 1, 1); book(_h_bjet_DY20 , 9, 1, 1); book(_h_bjet_ZdPhi20 ,11, 1, 1); book(_h_bjet_ZdR20 ,13, 1, 1); book(_h_bjet_ZPt ,15, 1, 1); book(_h_bjet_ZY ,17, 1, 1); book(_h_2bjet_dR ,21, 1, 1); book(_h_2bjet_Mbb ,23, 1, 1); book(_h_2bjet_ZPt ,25, 1, 1); book(_h_2bjet_ZY ,27, 1, 1); } //========================================================================================== /// Perform the per-event analysis void analyze(const Event& e) { - + //--------------------------- - const double weight = 1.0; + // -- check we have a Z: + const ZFinder& zfinder = apply(e, "ZFinder"); - // -- check we have a Z: - const ZFinder& zfinder = apply(e, "ZFinder"); - if(zfinder.bosons().size() != 1) vetoEvent; - + const ParticleVector boson_s = zfinder.bosons(); const Particle boson_f = boson_s[0] ; const ParticleVector zleps = zfinder.constituents(); //--------------------------- - + //--------------------------- //------------- stop processing the event if no true b-partons or hadrons are found const Particles& allBs = apply(e, "BHadrons").bHadrons(5.0*GeV); Particles stableBs; foreach(Particle p, allBs) { if(p.abseta() < 2.5) stableBs += p; } if( stableBs.empty() ) vetoEvent; - + //--------------------------- // -- get the b-jets: const Jets& jets = apply(e, "AntiKtJets04").jetsByPt(Cuts::pT >20.0*GeV && Cuts::abseta <2.4); Jets b_jets; foreach(const Jet& jet, jets) { //veto overlaps with Z leptons: bool veto = false; foreach(const Particle& zlep, zleps) { if(deltaR(jet, zlep) < 0.5) veto = true; } if(veto) continue; - + foreach(const Particle& bhadron, stableBs) { if( deltaR(jet, bhadron) <= 0.3 ) { b_jets.push_back(jet); break; // match } - } // end loop on b-hadrons + } // end loop on b-hadrons } - + //and make sure we have at least 1: if(b_jets.empty()) vetoEvent; //--------------------------- // fill the plots: const double ZpT = boson_f.pT()/GeV; const double ZY = boson_f.absrap(); - - _h_bjet_ZPt->fill(ZpT, weight); - _h_bjet_ZY ->fill(ZY, weight); - + + _h_bjet_ZPt->fill(ZpT); + _h_bjet_ZY ->fill(ZY); + foreach(const Jet& jet, b_jets) { - - _h_bjet_Pt->fill(jet.pT()/GeV, weight ); - _h_bjet_Y ->fill(jet.absrap(), weight ); - + + _h_bjet_Pt->fill(jet.pT()/GeV); + _h_bjet_Y ->fill(jet.absrap()); + const double Yboost = 0.5 * fabs(boson_f.rapidity() + jet.rapidity()); - _h_bjet_Yboost->fill(Yboost, weight ); - + _h_bjet_Yboost->fill(Yboost); + if(ZpT > 20.) { - + const double ZBDY = fabs( boson_f.rapidity() - jet.rapidity() ); const double ZBDPHI = fabs( deltaPhi(jet.phi(), boson_f.phi()) ); const double ZBDR = deltaR(jet, boson_f, RAPIDITY); - _h_bjet_DY20->fill( ZBDY, weight); - _h_bjet_ZdPhi20->fill(ZBDPHI, weight); - _h_bjet_ZdR20->fill( ZBDR, weight); + _h_bjet_DY20->fill( ZBDY); + _h_bjet_ZdPhi20->fill(ZBDPHI); + _h_bjet_ZdR20->fill( ZBDR); } - + } //loop over b-jets - + if (b_jets.size() < 2) return; - - _h_2bjet_ZPt->fill(ZpT, weight); - _h_2bjet_ZY ->fill(ZY, weight); - + + _h_2bjet_ZPt->fill(ZpT); + _h_2bjet_ZY ->fill(ZY); + const double BBDR = deltaR(b_jets[0], b_jets[1], RAPIDITY); const double Mbb = (b_jets[0].momentum() + b_jets[1].momentum()).mass(); - - _h_2bjet_dR ->fill(BBDR, weight); - _h_2bjet_Mbb->fill(Mbb, weight); - + + _h_2bjet_dR ->fill(BBDR); + _h_2bjet_Mbb->fill(Mbb); + } // end of analysis loop - - + + /// Normalise histograms etc., after the run void finalize() { - + const double normfac = crossSection() / sumOfWeights(); - + scale( _h_bjet_Pt, normfac); scale( _h_bjet_Y, normfac); scale( _h_bjet_Yboost, normfac); scale( _h_bjet_DY20, normfac); scale( _h_bjet_ZdPhi20, normfac); scale( _h_bjet_ZdR20, normfac); scale( _h_bjet_ZPt, normfac); scale( _h_bjet_ZY, normfac); scale( _h_2bjet_dR, normfac); scale( _h_2bjet_Mbb, normfac); scale( _h_2bjet_ZPt, normfac); scale( _h_2bjet_ZY, normfac); } - + //@} - - + + protected: - + // Data members like post-cuts event weight counters go here size_t _mode; - - + + private: Histo1DPtr _h_bjet_Pt; Histo1DPtr _h_bjet_Y; Histo1DPtr _h_bjet_Yboost; Histo1DPtr _h_bjet_DY20; Histo1DPtr _h_bjet_ZdPhi20; Histo1DPtr _h_bjet_ZdR20; Histo1DPtr _h_bjet_ZPt; Histo1DPtr _h_bjet_ZY; Histo1DPtr _h_2bjet_dR; Histo1DPtr _h_2bjet_Mbb; Histo1DPtr _h_2bjet_ZPt; Histo1DPtr _h_2bjet_ZY; - + }; class ATLAS_2014_I1306294_EL : public ATLAS_2014_I1306294 { public: ATLAS_2014_I1306294_EL() : ATLAS_2014_I1306294("ATLAS_2014_I1306294_EL") { _mode = 1; } }; class ATLAS_2014_I1306294_MU : public ATLAS_2014_I1306294 { public: ATLAS_2014_I1306294_MU() : ATLAS_2014_I1306294("ATLAS_2014_I1306294_MU") { _mode = 2; } }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2014_I1306294); DECLARE_RIVET_PLUGIN(ATLAS_2014_I1306294_MU); DECLARE_RIVET_PLUGIN(ATLAS_2014_I1306294_EL); -} - +} diff --git a/analyses/pluginATLAS/ATLAS_2014_I1306615.cc b/analyses/pluginATLAS/ATLAS_2014_I1306615.cc --- a/analyses/pluginATLAS/ATLAS_2014_I1306615.cc +++ b/analyses/pluginATLAS/ATLAS_2014_I1306615.cc @@ -1,486 +1,482 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { /// @brief ATLAS H->yy differential cross-sections measurement /// /// @author Michaela Queitsch-Maitland // // arXiv: http://arxiv.org/abs/ARXIV:1407.4222 // HepData: http://hepdata.cedar.ac.uk/view/ins1306615 class ATLAS_2014_I1306615 : public Analysis { public: // Constructor ATLAS_2014_I1306615() : Analysis("ATLAS_2014_I1306615") { } // Book histograms and initialise projections before the run void init() { // Final state // All particles within |eta| < 5.0 const FinalState FS(Cuts::abseta<5.0); declare(FS,"FS"); // Project photons with pT > 25 GeV and |eta| < 2.37 IdentifiedFinalState ph_FS(Cuts::abseta<2.37 && Cuts::pT>25.0*GeV); ph_FS.acceptIdPair(PID::PHOTON); declare(ph_FS, "PH_FS"); // Project photons for dressing IdentifiedFinalState ph_dressing_FS(FS); ph_dressing_FS.acceptIdPair(PID::PHOTON); // Project bare electrons IdentifiedFinalState el_bare_FS(FS); el_bare_FS.acceptIdPair(PID::ELECTRON); declare(el_bare_FS,"el_bare_FS"); // Project dressed electrons with pT > 15 GeV and |eta| < 2.47 DressedLeptons el_dressed_FS(ph_dressing_FS, el_bare_FS, 0.1, Cuts::abseta < 2.47 && Cuts::pT > 15*GeV); declare(el_dressed_FS,"EL_DRESSED_FS"); // Project bare muons IdentifiedFinalState mu_bare_FS(FS); mu_bare_FS.acceptIdPair(PID::MUON); // Project dressed muons with pT > 15 GeV and |eta| < 2.47 //DressedLeptons mu_dressed_FS(ph_dressing_FS, mu_bare_FS, 0.1, true, -2.47, 2.47, 15.0*GeV, false); DressedLeptons mu_dressed_FS(ph_dressing_FS, mu_bare_FS, 0.1, Cuts::abseta < 2.47 && Cuts::pT > 15*GeV); declare(mu_dressed_FS,"MU_DRESSED_FS"); // Final state excluding muons and neutrinos (for jet building and photon isolation) VetoedFinalState veto_mu_nu_FS(FS); veto_mu_nu_FS.vetoNeutrinos(); veto_mu_nu_FS.addVetoPairId(PID::MUON); declare(veto_mu_nu_FS, "VETO_MU_NU_FS"); // Build the anti-kT R=0.4 jets, using FinalState particles (vetoing muons and neutrinos) FastJets jets(veto_mu_nu_FS, FastJets::ANTIKT, 0.4); declare(jets, "JETS"); // Book histograms // 1D distributions book(_h_pT_yy ,1,1,1); book(_h_y_yy ,2,1,1); book(_h_Njets30 ,3,1,1); book(_h_Njets50 ,4,1,1); book(_h_pT_j1 ,5,1,1); book(_h_y_j1 ,6,1,1); book(_h_HT ,7,1,1); book(_h_pT_j2 ,8,1,1); book(_h_Dy_jj ,9,1,1); book(_h_Dphi_yy_jj ,10,1,1); book(_h_cosTS_CS ,11,1,1); book(_h_cosTS_CS_5bin ,12,1,1); book(_h_Dphi_jj ,13,1,1); book(_h_pTt_yy ,14,1,1); book(_h_Dy_yy ,15,1,1); book(_h_tau_jet ,16,1,1); book(_h_sum_tau_jet ,17,1,1); book(_h_y_j2 ,18,1,1); book(_h_pT_j3 ,19,1,1); book(_h_m_jj ,20,1,1); book(_h_pT_yy_jj ,21,1,1); // 2D distributions of cosTS_CS x pT_yy book(_h_cosTS_pTyy_low ,22,1,1); book(_h_cosTS_pTyy_high ,22,1,2); book(_h_cosTS_pTyy_rest ,22,1,3); // 2D distributions of Njets x pT_yy book(_h_pTyy_Njets0 ,23,1,1); book(_h_pTyy_Njets1 ,23,1,2); book(_h_pTyy_Njets2 ,23,1,3); book(_h_pTj1_excl ,24,1,1); // Fiducial regions book(_h_fidXSecs ,29,1,1); } // Perform the per-event analysis void analyze(const Event& event) { - const double weight = 1.0; - _weight = weight; - // Get final state particles const ParticleVector& FS_ptcls = apply(event, "FS").particles(); const ParticleVector& ptcls_veto_mu_nu = apply(event, "VETO_MU_NU_FS").particles(); const ParticleVector& photons = apply(event, "PH_FS").particlesByPt(); const vector& el_dressed = apply(event, "EL_DRESSED_FS").dressedLeptons(); const vector& mu_dressed = apply(event, "MU_DRESSED_FS").dressedLeptons(); // For isolation calculation float dR_iso = 0.4; float ETcut_iso = 14.0; FourMomentum ET_iso; // Fiducial selection: pT > 25 GeV, |eta| < 2.37 and isolation (in cone deltaR = 0.4) is < 14 GeV vector fid_photons; foreach (const Particle& ph, photons) { // Veto photons from hadron or tau decay if ( fromHadronDecay(ph) ) continue; // Calculate isolation ET_iso = - ph.momentum(); // Loop over fs truth particles (excluding muons and neutrinos) foreach (const Particle& p, ptcls_veto_mu_nu) { // Check if the truth particle is in a cone of 0.4 if ( deltaR(ph.momentum(), p.momentum()) < dR_iso ) ET_iso += p.momentum(); } // Check isolation if ( ET_iso.Et() > ETcut_iso ) continue; // Fill vector of photons passing fiducial selection fid_photons.push_back(&ph); } if(fid_photons.size() < 2) vetoEvent; const FourMomentum& y1 = fid_photons[0]->momentum(); const FourMomentum& y2 = fid_photons[1]->momentum(); double m_yy = (y1 + y2).mass(); // Relative pT cuts if ( y1.pT() < 0.35 * m_yy || y2.pT() < 0.25 * m_yy ) vetoEvent; // Mass window cut if ( m_yy < 105 || m_yy > 160 ) vetoEvent; // -------------------------------------------- // // Passed diphoton baseline fiducial selection! // // -------------------------------------------- // // Electron selection vector good_el; foreach(const DressedLepton& els, el_dressed) { const Particle& el = els.constituentLepton(); if ( el.momentum().pT() < 15 ) continue; if ( fabs(el.momentum().eta()) > 2.47 ) continue; if ( deltaR(el.momentum(), y1) < 0.4 ) continue; if ( deltaR(el.momentum(), y2) < 0.4 ) continue; if ( fromHadronDecay(el) ) continue; // Veto electrons from hadron or tau decay good_el.push_back(&el); } // Muon selection vector good_mu; foreach(const DressedLepton& mus, mu_dressed) { const Particle& mu = mus.constituentLepton(); if ( mu.momentum().pT() < 15 ) continue; if ( fabs(mu.momentum().eta()) > 2.47 ) continue; if ( deltaR(mu.momentum(), y1) < 0.4 ) continue; if ( deltaR(mu.momentum(), y2) < 0.4 ) continue; if ( fromHadronDecay(mu) ) continue; // Veto muons from hadron or tau decay good_mu.push_back(&mu); } // Find prompt, invisible particles for missing ET calculation // Based on VisibleFinalState projection FourMomentum invisible(0,0,0,0); foreach (const Particle& p, FS_ptcls) { // Veto non-prompt particles (from hadron or tau decay) if ( fromHadronDecay(p) ) continue; // Charged particles are visible if ( PID::threeCharge( p.pid() ) != 0 ) continue; // Neutral hadrons are visible if ( PID::isHadron( p.pid() ) ) continue; // Photons are visible if ( p.pid() == PID::PHOTON ) continue; // Gluons are visible (for parton level analyses) if ( p.pid() == PID::GLUON ) continue; // Everything else is invisible invisible += p.momentum(); } double MET = invisible.Et(); // Jet selection // Get jets with pT > 25 GeV and |rapidity| < 4.4 //const Jets& jets = apply(event, "JETS").jetsByPt(25.0*GeV, MAXDOUBLE, -4.4, 4.4, RAPIDITY); const Jets& jets = apply(event, "JETS").jetsByPt(Cuts::pT>25.0*GeV && Cuts::absrap <4.4); vector jets_25; vector jets_30; vector jets_50; foreach (const Jet& jet, jets) { bool passOverlap = true; // Overlap with leading photons if ( deltaR(y1, jet.momentum()) < 0.4 ) passOverlap = false; if ( deltaR(y2, jet.momentum()) < 0.4 ) passOverlap = false; // Overlap with good electrons foreach (const Particle* el, good_el) if ( deltaR(el->momentum(), jet.momentum()) < 0.2 ) passOverlap = false; if ( ! passOverlap ) continue; if ( fabs(jet.momentum().eta()) < 2.4 || ( fabs(jet.momentum().eta()) > 2.4 && jet.momentum().pT() > 30 ) ) jets_25.push_back(&jet); if ( jet.momentum().pT() > 30 ) jets_30.push_back(&jet); if ( jet.momentum().pT() > 50 ) jets_50.push_back(&jet); } // Fiducial regions - _h_fidXSecs->fill(1,_weight); - if ( jets_30.size() >= 1 ) _h_fidXSecs->fill(2, _weight); - if ( jets_30.size() >= 2 ) _h_fidXSecs->fill(3, _weight); - if ( jets_30.size() >= 3 ) _h_fidXSecs->fill(4, _weight); - if ( jets_30.size() >= 2 && passVBFCuts(y1 + y2, jets_30.at(0)->momentum(), jets_30.at(1)->momentum()) ) _h_fidXSecs->fill(5, _weight); - if ( (good_el.size() + good_mu.size()) > 0 ) _h_fidXSecs->fill(6, _weight); - if ( MET > 80 ) _h_fidXSecs->fill(7, _weight); + _h_fidXSecs->fill(1); + if ( jets_30.size() >= 1 ) _h_fidXSecs->fill(2); + if ( jets_30.size() >= 2 ) _h_fidXSecs->fill(3); + if ( jets_30.size() >= 3 ) _h_fidXSecs->fill(4); + if ( jets_30.size() >= 2 && passVBFCuts(y1 + y2, jets_30.at(0)->momentum(), jets_30.at(1)->momentum()) ) _h_fidXSecs->fill(5); + if ( (good_el.size() + good_mu.size()) > 0 ) _h_fidXSecs->fill(6); + if ( MET > 80 ) _h_fidXSecs->fill(7); // Fill histograms // Inclusive variables _pT_yy = (y1 + y2).pT(); _y_yy = fabs( (y1 + y2).rapidity() ); _cosTS_CS = cosTS_CS(y1, y2); _pTt_yy = pTt(y1, y2); _Dy_yy = fabs( deltaRap(y1, y2) ); _Njets30 = jets_30.size() > 3 ? 3 : jets_30.size(); _Njets50 = jets_50.size() > 3 ? 3 : jets_50.size(); - _h_Njets30->fill(_Njets30, _weight); - _h_Njets50->fill(_Njets50, _weight); + _h_Njets30->fill(_Njets30); + _h_Njets50->fill(_Njets50); _pT_j1 = jets_30.size() > 0 ? jets_30.at(0)->momentum().pT() : 0.; _pT_j2 = jets_30.size() > 1 ? jets_30.at(1)->momentum().pT() : 0.; _pT_j3 = jets_30.size() > 2 ? jets_30.at(2)->momentum().pT() : 0.; _HT = 0.0; foreach (const Jet* jet, jets_30) _HT += jet->momentum().pT(); _tau_jet = tau_jet_max(y1 + y2, jets_25); _sum_tau_jet = sum_tau_jet(y1 + y2, jets_25); - _h_pT_yy ->fill(_pT_yy ,_weight); - _h_y_yy ->fill(_y_yy ,_weight); - _h_pT_j1 ->fill(_pT_j1 ,_weight); - _h_cosTS_CS ->fill(_cosTS_CS ,_weight); - _h_cosTS_CS_5bin->fill(_cosTS_CS ,_weight); - _h_HT ->fill(_HT ,_weight); - _h_pTt_yy ->fill(_pTt_yy ,_weight); - _h_Dy_yy ->fill(_Dy_yy ,_weight); - _h_tau_jet ->fill(_tau_jet ,_weight); - _h_sum_tau_jet ->fill(_sum_tau_jet,_weight); + _h_pT_yy ->fill(_pT_yy); + _h_y_yy ->fill(_y_yy); + _h_pT_j1 ->fill(_pT_j1); + _h_cosTS_CS ->fill(_cosTS_CS); + _h_cosTS_CS_5bin->fill(_cosTS_CS); + _h_HT ->fill(_HT); + _h_pTt_yy ->fill(_pTt_yy); + _h_Dy_yy ->fill(_Dy_yy); + _h_tau_jet ->fill(_tau_jet); + _h_sum_tau_jet ->fill(_sum_tau_jet); // >=1 jet variables if ( jets_30.size() >= 1 ) { FourMomentum j1 = jets_30[0]->momentum(); _y_j1 = fabs( j1.rapidity() ); - _h_pT_j2->fill(_pT_j2 ,_weight); - _h_y_j1 ->fill(_y_j1 ,_weight); + _h_pT_j2->fill(_pT_j2); + _h_y_j1 ->fill(_y_j1); } // >=2 jet variables if ( jets_30.size() >= 2 ) { FourMomentum j1 = jets_30[0]->momentum(); FourMomentum j2 = jets_30[1]->momentum(); _Dy_jj = fabs( deltaRap(j1, j2) ); _Dphi_jj = fabs( deltaPhi(j1, j2) ); _Dphi_yy_jj = fabs( deltaPhi(y1 + y2, j1 + j2) ); _m_jj = (j1 + j2).mass(); _pT_yy_jj = (y1 + y2 + j1 + j2).pT(); _y_j2 = fabs( j2.rapidity() ); - _h_Dy_jj ->fill(_Dy_jj ,_weight); - _h_Dphi_jj ->fill(_Dphi_jj ,_weight); - _h_Dphi_yy_jj ->fill(_Dphi_yy_jj,_weight); - _h_m_jj ->fill(_m_jj ,_weight); - _h_pT_yy_jj ->fill(_pT_yy_jj ,_weight); - _h_pT_j3 ->fill(_pT_j3 ,_weight); - _h_y_j2 ->fill(_y_j2 ,_weight); + _h_Dy_jj ->fill(_Dy_jj); + _h_Dphi_jj ->fill(_Dphi_jj); + _h_Dphi_yy_jj ->fill(_Dphi_yy_jj); + _h_m_jj ->fill(_m_jj); + _h_pT_yy_jj ->fill(_pT_yy_jj); + _h_pT_j3 ->fill(_pT_j3); + _h_y_j2 ->fill(_y_j2); } // 2D distributions of cosTS_CS x pT_yy if ( _pT_yy < 80 ) - _h_cosTS_pTyy_low->fill(_cosTS_CS, _weight); + _h_cosTS_pTyy_low->fill(_cosTS_CS); else if ( _pT_yy > 80 && _pT_yy < 200 ) - _h_cosTS_pTyy_high->fill(_cosTS_CS,_weight); + _h_cosTS_pTyy_high->fill(_cosTS_CS); else if ( _pT_yy > 200 ) - _h_cosTS_pTyy_rest->fill(_cosTS_CS,_weight); + _h_cosTS_pTyy_rest->fill(_cosTS_CS); // 2D distributions of pT_yy x Njets if ( _Njets30 == 0 ) - _h_pTyy_Njets0->fill(_pT_yy, _weight); + _h_pTyy_Njets0->fill(_pT_yy); else if ( _Njets30 == 1 ) - _h_pTyy_Njets1->fill(_pT_yy, _weight); + _h_pTyy_Njets1->fill(_pT_yy); else if ( _Njets30 >= 2 ) - _h_pTyy_Njets2->fill(_pT_yy, _weight); + _h_pTyy_Njets2->fill(_pT_yy); - if ( _Njets30 == 1 ) _h_pTj1_excl->fill(_pT_j1, _weight); + if ( _Njets30 == 1 ) _h_pTj1_excl->fill(_pT_j1); } // Normalise histograms after the run void finalize() { const double xs = crossSectionPerEvent()/femtobarn; scale(_h_pT_yy, xs); scale(_h_y_yy, xs); scale(_h_pT_j1, xs); scale(_h_y_j1, xs); scale(_h_HT, xs); scale(_h_pT_j2, xs); scale(_h_Dy_jj, xs); scale(_h_Dphi_yy_jj, xs); scale(_h_cosTS_CS, xs); scale(_h_cosTS_CS_5bin, xs); scale(_h_Dphi_jj, xs); scale(_h_pTt_yy, xs); scale(_h_Dy_yy, xs); scale(_h_tau_jet, xs); scale(_h_sum_tau_jet, xs); scale(_h_y_j2, xs); scale(_h_pT_j3, xs); scale(_h_m_jj, xs); scale(_h_pT_yy_jj, xs); scale(_h_cosTS_pTyy_low, xs); scale(_h_cosTS_pTyy_high, xs); scale(_h_cosTS_pTyy_rest, xs); scale(_h_pTyy_Njets0, xs); scale(_h_pTyy_Njets1, xs); scale(_h_pTyy_Njets2, xs); scale(_h_pTj1_excl, xs); scale(_h_Njets30, xs); scale(_h_Njets50, xs); scale(_h_fidXSecs, xs); } // Trace event record to see if particle came from a hadron (or a tau from a hadron decay) // Based on fromDecay() function bool fromHadronDecay(const Particle& p ) { const GenVertex* prodVtx = p.genParticle()->production_vertex(); if (prodVtx == NULL) return false; foreach (const GenParticle* ancestor, particles(prodVtx, HepMC::ancestors)) { const PdgId pid = ancestor->pdg_id(); if (ancestor->status() == 2 && PID::isHadron(pid)) return true; if (ancestor->status() == 2 && (abs(pid) == PID::TAU && fromHadronDecay(ancestor))) return true; } return false; } // VBF-enhanced dijet topology selection cuts bool passVBFCuts(const FourMomentum &H, const FourMomentum &j1, const FourMomentum &j2) { return ( fabs(deltaRap(j1, j2)) > 2.8 && (j1 + j2).mass() > 400 && fabs(deltaPhi(H, j1 + j2)) > 2.6 ); } // Cosine of the decay angle in the Collins-Soper frame double cosTS_CS(const FourMomentum &y1, const FourMomentum &y2) { return fabs( ( (y1.E() + y1.pz())* (y2.E() - y2.pz()) - (y1.E() - y1.pz()) * (y2.E() + y2.pz()) ) / ((y1 + y2).mass() * sqrt(pow((y1 + y2).mass(), 2) + pow((y1 + y2).pt(), 2)) ) ); } // Diphoton pT along thrust axis double pTt(const FourMomentum &y1, const FourMomentum &y2) { return fabs(y1.px() * y2.py() - y2.px() * y1.py()) / (y1 - y2).pT()*2; } // Tau of jet (see paper for description) // tau_jet = mT/(2*cosh(y*)), where mT = pT (+) m, and y* = rapidty in Higgs rest frame double tau_jet( const FourMomentum &H, const FourMomentum &jet ) { return sqrt( pow(jet.pT(),2) + pow(jet.mass(),2) ) / (2.0 * cosh( jet.rapidity() - H.rapidity() ) ); } // Maximal (leading) tau_jet (see paper for description) double tau_jet_max( const FourMomentum &H, const vector jets, double tau_jet_cut = 8. ) { double max_tj = 0; for (size_t i=0; i < jets.size(); ++i) { FourMomentum jet = jets[i]->momentum(); if ( tau_jet(H, jet) > tau_jet_cut ) max_tj = max( tau_jet(H, jet), max_tj ); } return max_tj; } // Scalar sum of tau for all jets (see paper for description) double sum_tau_jet( const FourMomentum &H, const vector jets, double tau_jet_cut = 8. ) { double sum_tj = 0; for (size_t i=0; i < jets.size(); ++i) { FourMomentum jet = jets[i]->momentum(); if ( tau_jet(H, jet) > tau_jet_cut ) sum_tj += tau_jet(H, jet); } return sum_tj; } private: Histo1DPtr _h_pT_yy; Histo1DPtr _h_y_yy; Histo1DPtr _h_Njets30; Histo1DPtr _h_Njets50; Histo1DPtr _h_pT_j1; Histo1DPtr _h_y_j1; Histo1DPtr _h_HT; Histo1DPtr _h_pT_j2; Histo1DPtr _h_Dy_jj; Histo1DPtr _h_Dphi_yy_jj; Histo1DPtr _h_cosTS_CS; Histo1DPtr _h_cosTS_CS_5bin; Histo1DPtr _h_Dphi_jj; Histo1DPtr _h_pTt_yy; Histo1DPtr _h_Dy_yy; Histo1DPtr _h_tau_jet; Histo1DPtr _h_sum_tau_jet; Histo1DPtr _h_y_j2; Histo1DPtr _h_pT_j3; Histo1DPtr _h_m_jj; Histo1DPtr _h_pT_yy_jj; Histo1DPtr _h_cosTS_pTyy_low; Histo1DPtr _h_cosTS_pTyy_high; Histo1DPtr _h_cosTS_pTyy_rest; Histo1DPtr _h_pTyy_Njets0; Histo1DPtr _h_pTyy_Njets1; Histo1DPtr _h_pTyy_Njets2; Histo1DPtr _h_pTj1_excl; Histo1DPtr _h_fidXSecs; - double _weight; int _Njets30; int _Njets50; double _pT_yy; double _y_yy; double _cosTS_CS; double _pT_j1; double _m_jj; double _y_j1; double _HT; double _pT_j2; double _y_j2; double _Dphi_yy_jj; double _pT_yy_jj; double _Dphi_jj; double _Dy_jj; double _pT_j3; double _pTt_yy; double _Dy_yy; double _tau_jet; double _sum_tau_jet; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2014_I1306615); } diff --git a/analyses/pluginATLAS/ATLAS_2014_I1307756.cc b/analyses/pluginATLAS/ATLAS_2014_I1307756.cc --- a/analyses/pluginATLAS/ATLAS_2014_I1307756.cc +++ b/analyses/pluginATLAS/ATLAS_2014_I1307756.cc @@ -1,154 +1,154 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { class ATLAS_2014_I1307756 : public Analysis { public: /// Constructor ATLAS_2014_I1307756() : Analysis("ATLAS_2014_I1307756") { } /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { /// Initialise and register projections here FinalState fs; declare(fs, "FS"); FastJets fj(fs, FastJets::KT, 0.5); fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec())); declare(fj, "KtJetsD05"); IdentifiedFinalState photonfs(Cuts::abseta < 2.37 && Cuts::pT > 22*GeV); photonfs.acceptId(PID::PHOTON); declare(photonfs, "photons"); // Initialize event count here: - _fidWeights = 0.; + book(_fidWeights, "_fidWeights"); } int getEtaBin(double eta) const { double aeta = fabs(eta); return binIndex(aeta, _eta_bins_areaoffset); } /// Perform the per-event analysis void analyze(const Event& event) { /// Require at least 2 photons in final state Particles photons = apply(event, "photons").particlesByPt(); if (photons.size() < 2) vetoEvent; // Get jet pT densities vector< vector > ptDensities(_eta_bins_areaoffset.size()-1); const auto clust_seq_area = apply(event, "KtJetsD05").clusterSeqArea(); for (const Jet& jet : apply(event, "KtJetsD05").jets()) { const double area = clust_seq_area->area(jet); if (area > 1e-4 && jet.abseta() < _eta_bins_areaoffset.back()) { ptDensities.at(getEtaBin(jet.abseta())) += jet.pT()/area; } } /// Compute the median energy density per eta bin vector ptDensity; for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; ++b) { ptDensity += ptDensities[b].empty() ? 0 : median(ptDensities[b]); } // Loop over photons and find isolated ones Particles isolated_photons; for (const Particle& ph : photons) { Particles fs = apply(event, "FS").particles(); FourMomentum mom_in_EtCone; for (const Particle& p : fs) { // Reject if the particle is not in DR=0.4 cone if (deltaR(ph.momentum(), p.momentum()) > 0.4) continue; // Reject if the particle falls in the photon core if (fabs(ph.eta() - p.eta()) < 0.025 * 7 * 0.5 && fabs(ph.phi() - p.phi()) < PI/128. * 5 * 0.5) continue; // Reject if the particle is a neutrino (muons are kept) if (p.isNeutrino()) continue; // Sum momenta mom_in_EtCone += p.momentum(); } // Subtract the UE correction (area*density) const double ETCONE_AREA = M_PI*.4*.4 - (7.0*.025)*(5.0*M_PI/128.); const double correction = ptDensity[getEtaBin(ph.eta())] * ETCONE_AREA; // Add isolated photon to list if (mom_in_EtCone.Et() - correction > 12*GeV) continue; isolated_photons.push_back(ph); } // Require at least two isolated photons if (isolated_photons.size() < 2) vetoEvent ; // Select leading pT pair std::sort(isolated_photons.begin(), isolated_photons.end(), cmpMomByPt); const FourMomentum& y1 = isolated_photons[0].momentum(); const FourMomentum& y2 = isolated_photons[1].momentum(); // Compute invariant mass const FourMomentum yy = y1 + y2; const double Myy = yy.mass(); // If Myy >= 110 GeV, apply relative cuts if (Myy >= 110*GeV && (y1.Et()/Myy < 0.4 || y2.Et()/Myy < 0.3) ) vetoEvent; // Add to cross-section - _fidWeights += 1.0; + _fidWeights->fill(); } /// @todo Add to the YODA output rather than print to log void finalize() { // Compute selection efficiency & statistical error - const double eff = _fidWeights/sumOfWeights(); + const double eff = _fidWeights->val()/sumOfWeights(); const double err = sqrt(eff*(1-eff)/numEvents()); // Compute fiducial cross-section in fb const double fidCrossSection = eff * crossSection()/femtobarn; // Print out result MSG_INFO("=================================================="); MSG_INFO("==== Total cross-section: " << crossSection()/femtobarn<< " fb"); MSG_INFO("==== Fiducial cross-section: " << fidCrossSection << " fb"); MSG_INFO("=================================================="); MSG_INFO("==== Selection efficiency: " << eff << " +/- " << err << " (statistical error)"); MSG_INFO("=================================================="); } //@} private: const vector _eta_bins_areaoffset = {0.0, 1.5, 3.0}; - double _fidWeights; + CounterPtr _fidWeights; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2014_I1307756); } diff --git a/analyses/pluginATLAS/ATLAS_2014_I1312627.cc b/analyses/pluginATLAS/ATLAS_2014_I1312627.cc --- a/analyses/pluginATLAS/ATLAS_2014_I1312627.cc +++ b/analyses/pluginATLAS/ATLAS_2014_I1312627.cc @@ -1,258 +1,257 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/WFinder.hh" #include "Rivet/Projections/ZFinder.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/VetoedFinalState.hh" namespace Rivet { /// Measurement of V+jets distributions, taken as ratios between W and Z channels class ATLAS_2014_I1312627 : public Analysis { public: /// @name Plotting helper types //@{ struct Plots { string ref; Histo1DPtr comp[2]; // (de)nominator components Scatter2DPtr ratio; // Rjets plot }; typedef map PlotMap; typedef PlotMap::value_type PlotMapPair; //@} /// Constructor ATLAS_2014_I1312627(std::string name="ATLAS_2014_I1312627") : Analysis(name) { _mode = 0; // using electron channel for combined data by default setNeedsCrossSection(true); } /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Set up cuts Cut cuts; if (_mode == 2) { // muon channel cuts = (Cuts::pT > 25.0*GeV) & Cuts::etaIn(-2.4, 2.4); } else if (_mode) { // electron channel cuts = (Cuts::pT > 25.0*GeV) & ( Cuts::etaIn(-2.47, -1.52) | Cuts::etaIn(-1.37, 1.37) | Cuts::etaIn(1.52, 2.47) ); } else { // combined data extrapolated to common phase space cuts = (Cuts::pT > 25.0*GeV) & Cuts::etaIn(-2.5, 2.5); } // Boson finders FinalState fs; WFinder wfinder(fs, cuts, _mode > 1? PID::MUON : PID::ELECTRON, 40.0*GeV, MAXDOUBLE, 0.0*GeV, 0.1, WFinder::CLUSTERNODECAY, WFinder::NOTRACK, WFinder::TRANSMASS); declare(wfinder, "WF"); ZFinder zfinder(fs, cuts, _mode > 1? PID::MUON : PID::ELECTRON, 66.0*GeV, 116.0*GeV, 0.1, ZFinder::CLUSTERNODECAY, ZFinder::NOTRACK); declare(zfinder, "ZF"); // Jets VetoedFinalState jet_fs(fs); jet_fs.addVetoOnThisFinalState(getProjection("WF")); jet_fs.addVetoOnThisFinalState(getProjection("ZF")); FastJets jets(jet_fs, FastJets::ANTIKT, 0.4); jets.useInvisibles(true); declare(jets, "Jets"); // Book Rjets plots _suff = string("-y0") + to_str(_mode + 1); hInit("Njets_incl", "d01"); // inclusive number of jets hInit("Njets_excl", "d04"); // exclusive number of jets hInit("Pt1_N1incl", "d05"); // leading jet pT, at least 1 jet hInit("Pt1_N1excl", "d06"); // leading jet pT, exactly 1 jet hInit("Pt1_N2incl", "d07"); // leading jet pT, at least 2 jets hInit("Pt1_N3incl", "d08"); // leading jet pT, at least 3 jets hInit("Pt2_N2incl", "d09"); // subleading jet pT, at least 2 jets hInit("Pt3_N3incl", "d10"); // sub-subleading jet pT, at least 3 jets hInit("ST_N2incl", "d11"); // scalar jet pT sum, at least 2 jets hInit("ST_N2excl", "d12"); // scalar jet pT sum, exactly 2 jets hInit("ST_N3incl", "d13"); // scalar jet pT sum, at least 3 jets hInit("ST_N3excl", "d14"); // scalar jet pT sum, exactly 3 jets hInit("DR_N2incl", "d15"); // deltaR(j1, j2), at least 2 jets hInit("DPhi_N2incl", "d16"); // deltaPhi(j1, j2), at least 2 jets hInit("Mjj_N2incl", "d17"); // mjj, at least 2 jets hInit("Rap1_N1incl", "d18"); // leading jet rapidity, at least 1 jet hInit("Rap2_N2incl", "d19"); // subleading jet rapidity, at least 2 jets hInit("Rap3_N3incl", "d20"); // sub-subleading jet rapidity, at least 3 jets // Also book numerator and denominator for Rjets plots foreach (PlotMapPair& str_plot, _plots) { book(str_plot.second.comp[0] , str_plot.second.ref + "2" + _suff, *(str_plot.second.ratio) ); book(str_plot.second.comp[1] , str_plot.second.ref + "3" + _suff, *(str_plot.second.ratio) ); } } /// Perform the per-event analysis void analyze(const Event& event) { // Retrieve boson candidate const WFinder& wf = apply(event, "WF"); const ZFinder& zf = apply(event, "ZF"); if (wf.empty() && zf.empty()) vetoEvent; // Retrieve jets const JetAlg& jetfs = apply(event, "Jets"); Jets all_jets = jetfs.jetsByPt(Cuts::pT > 30*GeV && Cuts::absrap < 4.4); // Apply boson cuts and fill histograms - const double weight = 1.0; if (zf.size() == 2) { const Particles& leptons = zf.constituents(); if (oppSign(leptons[0], leptons[1]) && deltaR(leptons[0], leptons[1]) > 0.2) - fillPlots(leptons, all_jets, weight, 1); + fillPlots(leptons, all_jets, 1); } if (!wf.empty()) { const Particles& leptons = wf.constituentLeptons(); if (wf.constituentNeutrino().pT() > 25*GeV && wf.mT() > 40*GeV ) - fillPlots(leptons, all_jets, weight, 0); + fillPlots(leptons, all_jets, 0); } } /// Normalise histograms etc., after the run void finalize() { /// Normalise, scale and otherwise manipulate histograms here const double sf( crossSection() / sumOfWeights() ); foreach (PlotMapPair& str_plot, _plots) { scale(str_plot.second.comp[0], sf); scale(str_plot.second.comp[1], sf); divide(str_plot.second.comp[0], str_plot.second.comp[1], str_plot.second.ratio); } } //@} /// Analysis helper functions //@{ /// Histogram filling function, to avoid duplication - void fillPlots(const Particles& leptons, Jets& all_jets, const double& weight, int isZ) { + void fillPlots(const Particles& leptons, Jets& all_jets, int isZ) { // Do jet-lepton overlap removal Jets jets; foreach (const Jet& j, all_jets) { /// @todo A nice place for a lambda and any() logic when C++11 is available bool keep = true; foreach (const Particle& l, leptons) keep &= deltaR(j, l) > 0.5; if (keep) jets.push_back(j); } // Calculate jet ST const size_t njets = jets.size(); double ST = 0.0; // scalar pT sum of all selected jets for (size_t i = 0; i < njets; ++i) ST += jets[i].pT() * GeV; // Fill jet histos - _plots["Njets_excl"].comp[isZ]->fill(njets + 0.5, weight); + _plots["Njets_excl"].comp[isZ]->fill(njets + 0.5); for (size_t i = 0; i <= njets; ++i) - _plots["Njets_incl"].comp[isZ]->fill(i + 0.5, weight); + _plots["Njets_incl"].comp[isZ]->fill(i + 0.5); if (njets > 0) { const double pT1 = jets[0].pT()/GeV; const double rap1 = jets[0].absrap(); - _plots["Pt1_N1incl" ].comp[isZ]->fill(pT1, weight); - _plots["Rap1_N1incl"].comp[isZ]->fill(rap1, weight); + _plots["Pt1_N1incl" ].comp[isZ]->fill(pT1); + _plots["Rap1_N1incl"].comp[isZ]->fill(rap1); if (njets == 1) { - _plots["Pt1_N1excl"].comp[isZ]->fill(pT1, weight); + _plots["Pt1_N1excl"].comp[isZ]->fill(pT1); } else if (njets > 1) { const double pT2 = jets[1].pT()/GeV; const double rap2 = jets[1].absrap(); const double dR = deltaR(jets[0], jets[1]); const double dPhi = deltaPhi(jets[0], jets[1]); const double mjj = (jets[0].momentum() + jets[1].momentum()).mass()/GeV; - _plots["Pt1_N2incl" ].comp[isZ]->fill(pT1, weight); - _plots["Pt2_N2incl" ].comp[isZ]->fill(pT2, weight); - _plots["Rap2_N2incl"].comp[isZ]->fill(rap2, weight); - _plots["DR_N2incl" ].comp[isZ]->fill(dR, weight); - _plots["DPhi_N2incl"].comp[isZ]->fill(dPhi, weight); - _plots["Mjj_N2incl" ].comp[isZ]->fill(mjj, weight); - _plots["ST_N2incl" ].comp[isZ]->fill(ST, weight); + _plots["Pt1_N2incl" ].comp[isZ]->fill(pT1); + _plots["Pt2_N2incl" ].comp[isZ]->fill(pT2); + _plots["Rap2_N2incl"].comp[isZ]->fill(rap2); + _plots["DR_N2incl" ].comp[isZ]->fill(dR); + _plots["DPhi_N2incl"].comp[isZ]->fill(dPhi); + _plots["Mjj_N2incl" ].comp[isZ]->fill(mjj); + _plots["ST_N2incl" ].comp[isZ]->fill(ST); if (njets == 2) { - _plots["ST_N2excl"].comp[isZ]->fill(ST, weight); + _plots["ST_N2excl"].comp[isZ]->fill(ST); } else if (njets > 2) { const double pT3 = jets[2].pT()/GeV; const double rap3 = jets[2].absrap(); - _plots["Pt1_N3incl" ].comp[isZ]->fill(pT1, weight); - _plots["Pt3_N3incl" ].comp[isZ]->fill(pT3, weight); - _plots["Rap3_N3incl"].comp[isZ]->fill(rap3, weight); - _plots["ST_N3incl" ].comp[isZ]->fill(ST, weight); + _plots["Pt1_N3incl" ].comp[isZ]->fill(pT1); + _plots["Pt3_N3incl" ].comp[isZ]->fill(pT3); + _plots["Rap3_N3incl"].comp[isZ]->fill(rap3); + _plots["ST_N3incl" ].comp[isZ]->fill(ST); if (njets == 3) - _plots["ST_N3excl"].comp[isZ]->fill(ST, weight); + _plots["ST_N3excl"].comp[isZ]->fill(ST); } } } } /// Helper for histogram initialisation void hInit(string label, string ident) { string pre = ident + "-x0"; _plots[label].ref = pre; book(_plots[label].ratio, pre + "1" + _suff, true); } //@} protected: // Data members size_t _mode; string _suff; private: /// @name Map of histograms PlotMap _plots; }; /// Electron-specific version of the ATLAS_2014_I1312627 R-jets analysis class ATLAS_2014_I1312627_EL : public ATLAS_2014_I1312627 { public: ATLAS_2014_I1312627_EL() : ATLAS_2014_I1312627("ATLAS_2014_I1312627_EL") { _mode = 1; } }; /// Muon-specific version of the ATLAS_2014_I1312627 R-jets analysis class ATLAS_2014_I1312627_MU : public ATLAS_2014_I1312627 { public: ATLAS_2014_I1312627_MU() : ATLAS_2014_I1312627("ATLAS_2014_I1312627_MU") { _mode = 2; } }; // Hooks for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2014_I1312627); DECLARE_RIVET_PLUGIN(ATLAS_2014_I1312627_EL); DECLARE_RIVET_PLUGIN(ATLAS_2014_I1312627_MU); } diff --git a/analyses/pluginATLAS/ATLAS_2014_I1315949.cc b/analyses/pluginATLAS/ATLAS_2014_I1315949.cc --- a/analyses/pluginATLAS/ATLAS_2014_I1315949.cc +++ b/analyses/pluginATLAS/ATLAS_2014_I1315949.cc @@ -1,226 +1,225 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/ZFinder.hh" namespace Rivet { class ATLAS_2014_I1315949 : public Analysis { public: /// Constructor ATLAS_2014_I1315949() : Analysis("ATLAS_2014_I1315949") { } void init() { FinalState fs; ZFinder zfinder(fs, Cuts::abseta<2.4 && Cuts::pT>20.0*GeV, PID::MUON, 66*GeV, 116*GeV, 0.1, ZFinder::CLUSTERNODECAY); declare(zfinder, "ZFinder"); ChargedFinalState cfs( zfinder.remainingFinalState() ); declare(cfs, "cfs"); book(_h_pTsum_tow , 1, 1, 1); book(_h_pTsum_trv , 1, 1, 2); book(_h_pTsum_away , 1, 1, 3); book(_h_pTsum_tmin , 1, 1, 4); book(_h_pTsum_tmax , 1, 1, 5); book(_h_pTsum_tdif , 1, 1, 6); book(_h_Nchg_tow , 2, 1, 1); book(_h_Nchg_trv , 2, 1, 2); book(_h_Nchg_away , 2, 1, 3); book(_h_Nchg_tmin , 2, 1, 4); book(_h_Nchg_tmax , 2, 1, 5); book(_h_Nchg_tdif , 2, 1, 6); book(_h_pTavg_tow , 3, 1, 1); book(_h_pTavg_trv , 3, 1, 2); book(_h_pTavg_away , 3, 1, 3); book(_h_pTavgvsmult_tow , 4, 1, 1); book(_h_pTavgvsmult_trv , 4, 1, 2); book(_h_pTavgvsmult_away , 4, 1, 3); // Book sumpt and nch histos for (int i_reg = 0; i_reg < 4; i_reg++) { for (int i_bin = 0; i_bin < 6.; i_bin++) { book(_h_ptSum_1D[i_reg][i_bin], 5, i_reg+1 , i_bin+1); book( _h_Nchg_1D[i_reg][i_bin], 6, i_reg+1 , i_bin+1); } } } /// Perform the per-event analysis void analyze(const Event& event) { - const double weight = 1.0; const ZFinder& zfinder = apply(event, "ZFinder"); if (zfinder.bosons().size() != 1) vetoEvent; double Zpt = zfinder.bosons()[0].momentum().pT()/GeV; double Zphi = zfinder.bosons()[0].momentum().phi(); double Zmass = zfinder.bosons()[0].momentum().mass()/GeV; if(Zmass < 66. || Zmass > 116.) vetoEvent; // Initialise counters for Nch and sumPt for all regions int nTowards(0), nTransverse(0), nLeft(0), nRight(0), nTrmin(0), nTrmax(0), nAway(0); double ptSumTowards(0.0), ptSumTransverse(0.0), ptSumLeft(0.0), ptSumRight(0.0), ptSumTrmin(0.0), ptSumTrmax(0.0), ptSumAway(0.0); // The charged particles ParticleVector particles = apply(event, "cfs").particlesByPt( Cuts::pT > 0.5*GeV && Cuts::abseta <2.5); // Loop over charged particles with pT>500 MeV and |eta|<2.5 foreach(const Particle& p, particles) { double dphi = p.momentum().phi() - Zphi, pT = p.momentum().pT(); // Get multiples of 2pi right for(; std::fabs(dphi) > M_PI; dphi += (dphi > 0. ? -2.*M_PI : 2.*M_PI) ); // Towards region if( std::fabs(dphi) < M_PI/3. ) { nTowards++; ptSumTowards += pT; } // Transverse region else if( std::fabs(dphi) < 2.*M_PI/3. ) { nTransverse++; ptSumTransverse += pT; if(dphi > 0.) { nRight++; ptSumRight += pT; } else { nLeft++; ptSumLeft += pT; } } // Away region else { nAway++; ptSumAway += pT; } } // TransMAX, TransMIN regions if (ptSumLeft > ptSumRight) { ptSumTrmax = ptSumLeft; ptSumTrmin = ptSumRight; nTrmax = nLeft; nTrmin = nRight; } else { ptSumTrmax = ptSumRight; ptSumTrmin = ptSumLeft; nTrmax = nRight; nTrmin = nLeft; } // min max regions have difference are than all other regions const double area = 5.*2./3.*M_PI; // Fill sumPt vs. Zpt region profiles - _h_pTsum_tow->fill( Zpt, ptSumTowards/area, weight); - _h_pTsum_trv->fill( Zpt, ptSumTransverse/area, weight); - _h_pTsum_away->fill(Zpt, ptSumAway/area, weight); - _h_pTsum_tmin->fill(Zpt, ptSumTrmin/(0.5*area), weight); - _h_pTsum_tmax->fill(Zpt, ptSumTrmax/(0.5*area), weight); - _h_pTsum_tdif->fill(Zpt, (ptSumTrmax - ptSumTrmin)/(0.5*area), weight); + _h_pTsum_tow->fill( Zpt, ptSumTowards/area); + _h_pTsum_trv->fill( Zpt, ptSumTransverse/area); + _h_pTsum_away->fill(Zpt, ptSumAway/area); + _h_pTsum_tmin->fill(Zpt, ptSumTrmin/(0.5*area)); + _h_pTsum_tmax->fill(Zpt, ptSumTrmax/(0.5*area)); + _h_pTsum_tdif->fill(Zpt, (ptSumTrmax - ptSumTrmin)/(0.5*area)); // Fill Nch vs. Zpt region profiles - _h_Nchg_tow->fill( Zpt, nTowards/area, weight); - _h_Nchg_trv->fill( Zpt, nTransverse/area, weight); - _h_Nchg_away->fill(Zpt, nAway/area, weight); - _h_Nchg_tmin->fill(Zpt, nTrmin/(0.5*area), weight); - _h_Nchg_tmax->fill(Zpt, nTrmax/(0.5*area), weight); - _h_Nchg_tdif->fill(Zpt, (nTrmax - nTrmin)/(0.5*area), weight); + _h_Nchg_tow->fill( Zpt, nTowards/area); + _h_Nchg_trv->fill( Zpt, nTransverse/area); + _h_Nchg_away->fill(Zpt, nAway/area); + _h_Nchg_tmin->fill(Zpt, nTrmin/(0.5*area)); + _h_Nchg_tmax->fill(Zpt, nTrmax/(0.5*area)); + _h_Nchg_tdif->fill(Zpt, (nTrmax - nTrmin)/(0.5*area)); // Fill vs. ZpT profiles - _h_pTavg_tow->fill( Zpt, nTowards > 0.? ptSumTowards/nTowards : 0., weight); - _h_pTavg_trv->fill( Zpt, nTransverse > 0.? ptSumTransverse/nTransverse : 0., weight); - _h_pTavg_away->fill(Zpt, nAway > 0.? ptSumAway/nAway : 0., weight); + _h_pTavg_tow->fill( Zpt, nTowards > 0.? ptSumTowards/nTowards : 0.); + _h_pTavg_trv->fill( Zpt, nTransverse > 0.? ptSumTransverse/nTransverse : 0.); + _h_pTavg_away->fill(Zpt, nAway > 0.? ptSumAway/nAway : 0.); // Fill vs. ZpT profiles - _h_pTavgvsmult_tow->fill( nTowards, nTowards > 0.? ptSumTowards/nTowards : 0., weight); - _h_pTavgvsmult_trv->fill( nTransverse, nTransverse > 0.? ptSumTransverse/nTransverse : 0., weight); - _h_pTavgvsmult_away->fill(nAway, nAway > 0.? ptSumAway/nAway : 0., weight); + _h_pTavgvsmult_tow->fill( nTowards, nTowards > 0.? ptSumTowards/nTowards : 0.); + _h_pTavgvsmult_trv->fill( nTransverse, nTransverse > 0.? ptSumTransverse/nTransverse : 0.); + _h_pTavgvsmult_away->fill(nAway, nAway > 0.? ptSumAway/nAway : 0.); // Determine Zpt region histo to fill int i_bin(0); if (inRange(Zpt,0,5) ) i_bin=0; if (inRange(Zpt,5,10) ) i_bin=1; if (inRange(Zpt,10,20) ) i_bin=2; if (inRange(Zpt,20,50) ) i_bin=3; if (inRange(Zpt,50,110) ) i_bin=4; if (Zpt>110) i_bin=5; // SumPt histos for Zpt region - _h_ptSum_1D[0][i_bin]->fill(ptSumTowards/area, weight); - _h_ptSum_1D[1][i_bin]->fill(ptSumTransverse/area, weight); - _h_ptSum_1D[2][i_bin]->fill(ptSumTrmin/(0.5*area), weight); - _h_ptSum_1D[3][i_bin]->fill(ptSumTrmax/(0.5*area), weight); + _h_ptSum_1D[0][i_bin]->fill(ptSumTowards/area); + _h_ptSum_1D[1][i_bin]->fill(ptSumTransverse/area); + _h_ptSum_1D[2][i_bin]->fill(ptSumTrmin/(0.5*area)); + _h_ptSum_1D[3][i_bin]->fill(ptSumTrmax/(0.5*area)); // Nch histos for Zpt region - _h_Nchg_1D[0][i_bin]->fill(nTowards/area, weight); - _h_Nchg_1D[1][i_bin]->fill(nTransverse/area, weight); - _h_Nchg_1D[2][i_bin]->fill(nTrmin/(0.5*area), weight); - _h_Nchg_1D[3][i_bin]->fill(nTrmax/(0.5*area), weight); + _h_Nchg_1D[0][i_bin]->fill(nTowards/area); + _h_Nchg_1D[1][i_bin]->fill(nTransverse/area); + _h_Nchg_1D[2][i_bin]->fill(nTrmin/(0.5*area)); + _h_Nchg_1D[3][i_bin]->fill(nTrmax/(0.5*area)); } /// Normalise histograms etc., after the run void finalize() { for(int i_reg = 0; i_reg < 4; i_reg++) { for(int i_bin = 0; i_bin < 6; i_bin++) { normalize( _h_ptSum_1D[i_reg][i_bin] ); normalize( _h_Nchg_1D[ i_reg][i_bin] ); } } } private: Profile1DPtr _h_pTsum_tow, _h_pTsum_trv, _h_pTsum_away, _h_pTsum_tmin, _h_pTsum_tmax, _h_pTsum_tdif, _h_Nchg_tow, _h_Nchg_trv, _h_Nchg_away, _h_Nchg_tmin, _h_Nchg_tmax, _h_Nchg_tdif, _h_pTavg_tow, _h_pTavg_trv, _h_pTavg_away, _h_pTavgvsmult_tow, _h_pTavgvsmult_trv, _h_pTavgvsmult_away; Histo1DPtr _h_ptSum_1D[4][6], _h_Nchg_1D[4][6]; }; // This global object acts as a hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2014_I1315949); } diff --git a/analyses/pluginATLAS/ATLAS_2014_I1319490.cc b/analyses/pluginATLAS/ATLAS_2014_I1319490.cc --- a/analyses/pluginATLAS/ATLAS_2014_I1319490.cc +++ b/analyses/pluginATLAS/ATLAS_2014_I1319490.cc @@ -1,232 +1,232 @@ #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/WFinder.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/VetoedFinalState.hh" namespace Rivet { class ATLAS_2014_I1319490 : public Analysis { public: ATLAS_2014_I1319490(string name = "ATLAS_2014_I1319490") : Analysis(name) { _mode = 0; // using electron channel for combined data by default setNeedsCrossSection(true); } // Book histograms and initialise projections before the run void init() { FinalState fs; Cut cuts; if (_mode == 2) { // muon channel cuts = (Cuts::pT > 25.0*GeV) & Cuts::etaIn(-2.4, 2.4); } else if (_mode) { // electron channel cuts = (Cuts::pT > 25.0*GeV) & ( Cuts::etaIn(-2.47, -1.52) | Cuts::etaIn(-1.37, 1.37) | Cuts::etaIn(1.52, 2.47) ); } else { // combined data extrapolated to common phase space cuts = (Cuts::pT > 25.0*GeV) & Cuts::etaIn(-2.5, 2.5); } // bosons WFinder wfinder(fs, cuts, _mode > 1? PID::MUON : PID::ELECTRON, 40.0*GeV, MAXDOUBLE, 0.0*GeV, 0.1, WFinder::CLUSTERNODECAY, WFinder::NOTRACK, WFinder::TRANSMASS); declare(wfinder, "WF"); // jets VetoedFinalState jet_fs(fs); jet_fs.addVetoOnThisFinalState(getProjection("WF")); FastJets jets(jet_fs, FastJets::ANTIKT, 0.4); jets.useInvisibles(true); declare(jets, "Jets"); // book histograms book(histos["h_N_incl"] ,1,1,_mode+1); book(histos["h_N"] ,4,1,_mode+1); book(histos["h_pt_jet1_1jet"] ,5,1,_mode+1); book(histos["h_pt_jet1_1jet_excl"] ,6,1,_mode+1); book(histos["h_pt_jet1_2jet"] ,7,1,_mode+1); book(histos["h_pt_jet1_3jet"] ,8,1,_mode+1); book(histos["h_pt_jet2_2jet"] ,9,1,_mode+1); book(histos["h_pt_jet3_3jet"] ,10,1,_mode+1); book(histos["h_pt_jet4_4jet"] ,11,1,_mode+1); book(histos["h_pt_jet5_5jet"] ,12,1,_mode+1); book(histos["h_y_jet1_1jet"] ,13,1,_mode+1); book(histos["h_y_jet2_2jet"] ,14,1,_mode+1); book(histos["h_HT_1jet"] ,15,1,_mode+1); book(histos["h_HT_1jet_excl"] ,16,1,_mode+1); book(histos["h_HT_2jet"] ,17,1,_mode+1); book(histos["h_HT_2jet_excl"] ,18,1,_mode+1); book(histos["h_HT_3jet"] ,19,1,_mode+1); book(histos["h_HT_3jet_excl"] ,20,1,_mode+1); book(histos["h_HT_4jet"] ,21,1,_mode+1); book(histos["h_HT_5jet"] ,22,1,_mode+1); book(histos["h_deltaPhi_jet12"] ,23,1,_mode+1); book(histos["h_deltaRap_jet12"] ,24,1,_mode+1); book(histos["h_deltaR_jet12"] ,25,1,_mode+1); book(histos["h_M_Jet12_2jet"] ,26,1,_mode+1); book(histos["h_y_jet3_3jet"] ,27,1,_mode+1); book(histos["h_y_jet4_4jet"] ,28,1,_mode+1); book(histos["h_y_jet5_5jet"] ,29,1,_mode+1); book(histos["h_ST_1jet"] ,30,1,_mode+1); book(histos["h_ST_2jet"] ,31,1,_mode+1); book(histos["h_ST_2jet_excl"] ,32,1,_mode+1); book(histos["h_ST_3jet"] ,33,1,_mode+1); book(histos["h_ST_3jet_excl"] ,34,1,_mode+1); book(histos["h_ST_4jet"] ,35,1,_mode+1); book(histos["h_ST_5jet"] ,36,1,_mode+1); } - void fillPlots(const Particle& lepton, const double& missET, Jets& all_jets, const double& weight) { + void fillPlots(const Particle& lepton, const double& missET, Jets& all_jets) { // do jet-lepton overlap removal Jets jets; double ST = 0.0; // scalar pT sum of all selected jets foreach (const Jet &j, all_jets) { if (deltaR(j, lepton) > 0.5) { jets += j; ST += j.pT() / GeV; } } const size_t njets = jets.size(); const double HT = ST + lepton.pT() / GeV + missET; - histos["h_N"]->fill(njets + 0.5, weight); + histos["h_N"]->fill(njets + 0.5); for (size_t i = 0; i <= njets; ++i) { - histos["h_N_incl"]->fill(i + 0.5, weight); + histos["h_N_incl"]->fill(i + 0.5); } if (njets) { const double pT1 = jets[0].pT() / GeV; const double rap1 = jets[0].absrap(); - histos["h_pt_jet1_1jet" ]->fill(pT1, weight); - histos["h_y_jet1_1jet"]->fill(rap1, weight); - histos["h_HT_1jet"]->fill(HT, weight); - histos["h_ST_1jet"]->fill(ST, weight); + histos["h_pt_jet1_1jet" ]->fill(pT1); + histos["h_y_jet1_1jet"]->fill(rap1); + histos["h_HT_1jet"]->fill(HT); + histos["h_ST_1jet"]->fill(ST); if (njets == 1) { - histos["h_pt_jet1_1jet_excl"]->fill(pT1, weight); - histos["h_HT_1jet_excl"]->fill(HT, weight); + histos["h_pt_jet1_1jet_excl"]->fill(pT1); + histos["h_HT_1jet_excl"]->fill(HT); } else { const double pT2 = jets[1].pT() / GeV; const double rap2 = jets[1].absrap(); const double dR = deltaR(jets[0], jets[1]); const double dRap = deltaRap(jets[0], jets[1]); const double dPhi = deltaPhi(jets[0], jets[1]); const double mjj = (jets[0].momentum() + jets[1].momentum()).mass() / GeV; - histos["h_pt_jet1_2jet"]->fill(pT1, weight); - histos["h_pt_jet2_2jet"]->fill(pT2, weight); - histos["h_y_jet2_2jet"]->fill(rap2, weight); - histos["h_M_Jet12_2jet"]->fill(mjj, weight); - histos["h_HT_2jet"]->fill(HT, weight); - histos["h_ST_2jet"]->fill(ST, weight); - histos["h_deltaPhi_jet12"]->fill(dPhi, weight); - histos["h_deltaRap_jet12"]->fill(dRap, weight); - histos["h_deltaR_jet12"]->fill(dR, weight); + histos["h_pt_jet1_2jet"]->fill(pT1); + histos["h_pt_jet2_2jet"]->fill(pT2); + histos["h_y_jet2_2jet"]->fill(rap2); + histos["h_M_Jet12_2jet"]->fill(mjj); + histos["h_HT_2jet"]->fill(HT); + histos["h_ST_2jet"]->fill(ST); + histos["h_deltaPhi_jet12"]->fill(dPhi); + histos["h_deltaRap_jet12"]->fill(dRap); + histos["h_deltaR_jet12"]->fill(dR); if (njets == 2) { - histos["h_ST_2jet_excl"]->fill(ST, weight); - histos["h_HT_2jet_excl"]->fill(HT, weight); + histos["h_ST_2jet_excl"]->fill(ST); + histos["h_HT_2jet_excl"]->fill(HT); } else { const double pT3 = jets[2].pT() / GeV; const double rap3 = jets[2].absrap(); - histos["h_pt_jet1_3jet"]->fill(pT1, weight); - histos["h_pt_jet3_3jet"]->fill(pT3, weight); - histos["h_y_jet3_3jet"]->fill(rap3, weight); - histos["h_HT_3jet"]->fill(HT, weight); - histos["h_ST_3jet"]->fill(ST, weight); + histos["h_pt_jet1_3jet"]->fill(pT1); + histos["h_pt_jet3_3jet"]->fill(pT3); + histos["h_y_jet3_3jet"]->fill(rap3); + histos["h_HT_3jet"]->fill(HT); + histos["h_ST_3jet"]->fill(ST); if(njets == 3) { - histos["h_ST_3jet_excl"]->fill(ST, weight); - histos["h_HT_3jet_excl"]->fill(HT, weight); + histos["h_ST_3jet_excl"]->fill(ST); + histos["h_HT_3jet_excl"]->fill(HT); } else { const double pT4 = jets[3].pT() / GeV; const double rap4 = jets[3].absrap(); - histos["h_pt_jet4_4jet"]->fill(pT4, weight); - histos["h_y_jet4_4jet"]->fill(rap4, weight); - histos["h_HT_4jet"]->fill(HT, weight); - histos["h_ST_4jet"]->fill(ST, weight); + histos["h_pt_jet4_4jet"]->fill(pT4); + histos["h_y_jet4_4jet"]->fill(rap4); + histos["h_HT_4jet"]->fill(HT); + histos["h_ST_4jet"]->fill(ST); if (njets > 4) { const double pT5 = jets[4].pT() / GeV; const double rap5 = jets[4].absrap(); - histos["h_pt_jet5_5jet"]->fill(pT5, weight); - histos["h_y_jet5_5jet"]->fill(rap5, weight); - histos["h_HT_5jet"]->fill(HT, weight); - histos["h_ST_5jet"]->fill(ST, weight); + histos["h_pt_jet5_5jet"]->fill(pT5); + histos["h_y_jet5_5jet"]->fill(rap5); + histos["h_HT_5jet"]->fill(HT); + histos["h_ST_5jet"]->fill(ST); } } } } } } // Perform the per-event analysis void analyze(const Event& event) { // Retrieve boson candidate const WFinder& wf = apply(event, "WF"); if (wf.empty()) vetoEvent; // Retrieve jets const JetAlg& jetfs = apply(event, "Jets"); Jets all_jets = jetfs.jetsByPt(Cuts::pT > 30.0*GeV && Cuts::absrap < 4.4); const Particles& leptons = wf.constituentLeptons(); const double missET = wf.constituentNeutrino().pT() / GeV; if (leptons.size() == 1 && missET > 25.0 && wf.mT() > 40.0*GeV) { const Particle& lep = leptons[0]; - fillPlots(lep, missET, all_jets, 1.0); + fillPlots(lep, missET, all_jets); } } void finalize() { const double scalefactor(crossSection() / sumOfWeights()); /// @todo Update to use C++11 range-for for (map::iterator hit = histos.begin(); hit != histos.end(); ++hit) { scale(hit->second, scalefactor); } } protected: size_t _mode; private: map histos; }; class ATLAS_2014_I1319490_EL : public ATLAS_2014_I1319490 { public: ATLAS_2014_I1319490_EL() : ATLAS_2014_I1319490("ATLAS_2014_I1319490_EL") { _mode = 1; } }; class ATLAS_2014_I1319490_MU : public ATLAS_2014_I1319490 { public: ATLAS_2014_I1319490_MU() : ATLAS_2014_I1319490("ATLAS_2014_I1319490_MU") { _mode = 2; } }; // The hooks for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2014_I1319490); DECLARE_RIVET_PLUGIN(ATLAS_2014_I1319490_EL); DECLARE_RIVET_PLUGIN(ATLAS_2014_I1319490_MU); } diff --git a/analyses/pluginATLAS/ATLAS_2014_I1327229.cc b/analyses/pluginATLAS/ATLAS_2014_I1327229.cc --- a/analyses/pluginATLAS/ATLAS_2014_I1327229.cc +++ b/analyses/pluginATLAS/ATLAS_2014_I1327229.cc @@ -1,1330 +1,1326 @@ // -*- C++ -*- #include "Rivet/Analysis.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_2014_I1327229 : public Analysis { public: /// Constructor ATLAS_2014_I1327229() : Analysis("ATLAS_2014_I1327229") { } /// 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 /// @todo Replace with SmearedParticles etc. srand(160385); // Read in all signal regions _signal_regions = getSignalRegions(); // Set number of events per signal region to 0 for (size_t i = 0; i < _signal_regions.size(); i++) - _eventCountsPerSR[_signal_regions[i]] = 0.0; + book(_eventCountsPerSR[_signal_regions[i]], "_eventCountsPerSR_" + _signal_regions[i]); // Final state including all charged and neutral particles const FinalState fs(-5.0, 5.0, 1*GeV); declare(fs, "FS"); // Final state including all charged particles declare(ChargedFinalState(-2.5, 2.5, 1*GeV), "CFS"); // Final state including all visible particles (to calculate MET, Jets etc.) declare(VisibleFinalState(-5.0,5.0),"VFS"); // Final state including all AntiKt 04 Jets VetoedFinalState vfs; vfs.addVetoPairId(PID::MUON); declare(FastJets(vfs, FastJets::ANTIKT, 0.4), "AntiKtJets04"); // Final state including all unstable particles (including taus) declare(UnstableFinalState(Cuts::abseta < 5.0 && Cuts::pT > 5*GeV),"UFS"); // Final state including all electrons IdentifiedFinalState elecs(Cuts::abseta < 2.47 && Cuts::pT > 10*GeV); elecs.acceptIdPair(PID::ELECTRON); declare(elecs, "elecs"); // Final state including all muons IdentifiedFinalState muons(Cuts::abseta < 2.5 && Cuts::pT > 10*GeV); muons.acceptIdPair(PID::MUON); declare(muons, "muons"); /// Book histograms: book(_h_HTlep_all ,"HTlep_all", 30,0,3000); book(_h_HTjets_all ,"HTjets_all", 30,0,3000); book(_h_MET_all ,"MET_all", 30,0,1500); book(_h_Meff_all ,"Meff_all", 50,0,5000); book(_h_min_pT_all ,"min_pT_all", 50, 0, 2000); book(_h_mT_all ,"mT_all", 50, 0, 2000); book(_h_e_n ,"e_n", 10, -0.5, 9.5); book(_h_mu_n ,"mu_n", 10, -0.5, 9.5); book(_h_tau_n ,"tau_n", 10, -0.5, 9.5); book(_h_pt_1_3l ,"pt_1_3l", 100, 0, 2000); book(_h_pt_2_3l ,"pt_2_3l", 100, 0, 2000); book(_h_pt_3_3l ,"pt_3_3l", 100, 0, 2000); book(_h_pt_1_2ltau ,"pt_1_2ltau", 100, 0, 2000); book(_h_pt_2_2ltau ,"pt_2_2ltau", 100, 0, 2000); book(_h_pt_3_2ltau ,"pt_3_2ltau", 100, 0, 2000); book(_h_excluded ,"excluded", 2, -0.5, 1.5); } /// Perform the per-event analysis void analyze(const Event& event) { // Muons Particles muon_candidates; const Particles charged_tracks = apply(event, "CFS").particles(); const Particles visible_particles = apply(event, "VFS").particles(); for (const Particle& mu : apply(event, "muons").particlesByPt() ) { // Calculate pTCone30 variable (pT of all tracks within dR<0.3 - pT of muon itself) double pTinCone = -mu.pT(); for (const Particle& track : charged_tracks ) { if (deltaR(mu.momentum(),track.momentum()) < 0.3 ) pTinCone += track.pT(); } // Calculate eTCone30 variable (pT of all visible particles within dR<0.3) double eTinCone = 0.; for (const Particle& visible_particle : visible_particles) { if (visible_particle.abspid() != PID::MUON && inRange(deltaR(mu.momentum(),visible_particle.momentum()), 0.1, 0.3)) eTinCone += visible_particle.pT(); } // Apply reconstruction efficiency and simulate reconstruction int muon_id = 13; if (mu.hasAncestor(PID::TAU) || mu.hasAncestor(-PID::TAU)) muon_id = 14; const double eff = (_use_fiducial_lepton_efficiency) ? apply_reco_eff(muon_id,mu) : 1.0; const bool keep_muon = rand()/static_cast(RAND_MAX)<=eff; // Keep muon if pTCone30/pT < 0.15 and eTCone30/pT < 0.2 and reconstructed if (keep_muon && pTinCone/mu.pT() <= 0.1 && eTinCone/mu.pT() < 0.1) muon_candidates.push_back(mu); } // Electrons Particles electron_candidates; for (const Particle& e : apply(event, "elecs").particlesByPt() ) { // Neglect electrons in crack regions if (inRange(e.abseta(), 1.37, 1.52)) continue; // Calculate pTCone30 variable (pT of all tracks within dR<0.3 - pT of electron itself) double pTinCone = -e.pT(); for (const Particle& track : charged_tracks) { if (deltaR(e.momentum(), track.momentum()) < 0.3 ) pTinCone += track.pT(); } // Calculate eTCone30 variable (pT of all visible particles (except muons) within dR<0.3) double eTinCone = 0.; for (const Particle& visible_particle : visible_particles) { if (visible_particle.abspid() != PID::MUON && inRange(deltaR(e.momentum(),visible_particle.momentum()), 0.1, 0.3)) eTinCone += visible_particle.pT(); } // Apply reconstruction efficiency and simulate reconstruction int elec_id = 11; if (e.hasAncestor(15) || e.hasAncestor(-15)) elec_id = 12; const double eff = (_use_fiducial_lepton_efficiency) ? apply_reco_eff(elec_id,e) : 1.0; const bool keep_elec = rand()/static_cast(RAND_MAX)<=eff; // Keep electron if pTCone30/pT < 0.13 and eTCone30/pT < 0.2 and reconstructed if (keep_elec && pTinCone/e.pT() <= 0.1 && eTinCone/e.pT() < 0.1) electron_candidates.push_back(e); } // Taus Particles tau_candidates; for (const Particle& tau : apply(event, "UFS").particles() ) { // Only pick taus out of all unstable particles if ( tau.abspid() != PID::TAU) continue; // Check that tau has decayed into daughter particles if (tau.genParticle()->end_vertex() == 0) continue; // Calculate visible tau momentum using the tau neutrino momentum in the tau decay FourMomentum daughter_tau_neutrino_momentum = get_tau_neutrino_momentum(tau); Particle tau_vis = tau; tau_vis.setMomentum(tau.momentum()-daughter_tau_neutrino_momentum); // keep only taus in certain eta region and above 15 GeV of visible tau pT if ( tau_vis.pT()/GeV <= 15.0 || tau_vis.abseta() > 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); // Apply reconstruction efficiency and simulate reconstruction int tau_id = 15; if (nprong == 1) tau_id = 15; else if (nprong == 3) tau_id = 16; const double eff = (_use_fiducial_lepton_efficiency) ? apply_reco_eff(tau_id,tau_vis) : 1.0; const bool keep_tau = rand()/static_cast(RAND_MAX)<=eff; // Keep tau if nprong = 1, it decays hadronically and it is reconstructed if ( !lep_decaying_tau && nprong == 1 && keep_tau) tau_candidates.push_back(tau_vis); } // Jets (all anti-kt R=0.4 jets with pT > 30 GeV and eta < 4.9 Jets jet_candidates; for (const Jet& jet : apply(event, "AntiKtJets04").jetsByPt(30.0*GeV) ) { if (jet.abseta() < 4.9 ) jet_candidates.push_back(jet); } // ETmiss Particles vfs_particles = apply(event, "VFS").particles(); FourMomentum pTmiss; for (const Particle& p : vfs_particles) pTmiss -= p.momentum(); double eTmiss = pTmiss.pT()/GeV; // ------------------------- // Overlap removal // electron - electron Particles electron_candidates_2; for(size_t 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(size_t 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; for (const Jet& jet : jet_candidates) { bool away = true; // If jet within dR < 0.2 of electron: remove jet for (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 for (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 Particles recon_leptons, recon_e; for (size_t 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; for (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) { for (const Particle& mu : muon_candidates) { 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 Particles recon_tau; for (const Particle& tau : tau_candidates) { bool away = true; // If tau within dR < 0.2 of an electron: remove tau for (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) { for (const Particle& mu : muon_candidates) { if (deltaR(tau.momentum(), mu.momentum()) < 0.2 ) { away = false; break; } } } // If isolated keep it if (away) recon_tau.push_back( tau ); } // muon - jet Particles recon_mu, trigger_mu; // If muon within dR < 0.4 of a jet: remove muon for (const Particle& mu : muon_candidates ) { bool away = true; for (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 (mu.abseta() < 2.4) trigger_mu.push_back( mu ); } } // End overlap removal // --------------------- // Jet cleaning if (rand()/static_cast(RAND_MAX) <= 0.42) { for (const Jet& jet : recon_jets ) { const double eta = jet.rapidity(); const double phi = jet.azimuthalAngle(MINUSPI_PLUSPI); if(jet.pT() > 25*GeV && inRange(eta,-0.1,1.5) && inRange(phi,-0.9,-0.5)) vetoEvent; } } // Event selection // Require at least 3 charged tracks in event if (charged_tracks.size() < 3) vetoEvent; // And at least one e/mu passing trigger if( !( !recon_e.empty() && recon_e[0].pT()>26.*GeV) && !( !trigger_mu.empty() && trigger_mu[0].pT()>26.*GeV) ) { MSG_DEBUG("Hardest lepton fails trigger"); vetoEvent; } // And only accept 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; - // Getting the event weight - const double weight = 1.0; - // Sort leptons by decreasing pT sortByPt(recon_leptons); sortByPt(recon_tau); // Calculate HTlep, fill lepton pT histograms & store chosen combination of 3 leptons double HTlep = 0.; Particles chosen_leptons; if (recon_leptons.size() > 2) { - _h_pt_1_3l->fill(recon_leptons[0].pT()/GeV, weight); - _h_pt_2_3l->fill(recon_leptons[1].pT()/GeV, weight); - _h_pt_3_3l->fill(recon_leptons[2].pT()/GeV, weight); + _h_pt_1_3l->fill(recon_leptons[0].pT()/GeV); + _h_pt_2_3l->fill(recon_leptons[1].pT()/GeV); + _h_pt_3_3l->fill(recon_leptons[2].pT()/GeV); HTlep = (recon_leptons[0].pT() + recon_leptons[1].pT() + recon_leptons[2].pT())/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].pT()/GeV, weight); - _h_pt_2_2ltau->fill(recon_leptons[1].pT()/GeV, weight); - _h_pt_3_2ltau->fill(recon_tau[0].pT()/GeV, weight); + _h_pt_1_2ltau->fill(recon_leptons[0].pT()/GeV); + _h_pt_2_2ltau->fill(recon_leptons[1].pT()/GeV); + _h_pt_3_2ltau->fill(recon_tau[0].pT()/GeV); HTlep = recon_leptons[0].pT()/GeV + recon_leptons[1].pT()/GeV + recon_tau[0].pT()/GeV; chosen_leptons.push_back( recon_leptons[0] ); chosen_leptons.push_back( recon_leptons[1] ); chosen_leptons.push_back( recon_tau[0] ); } // Calculate mT and mTW variable Particles mT_leptons; Particles mTW_leptons; for (size_t i1 = 0; i1 < 3; i1 ++) { for (size_t i2 = i1+1; i2 < 3; i2 ++) { double OSSF_inv_mass = isOSSF_mass(chosen_leptons[i1],chosen_leptons[i2]); if (OSSF_inv_mass != 0.) { for (size_t i3 = 0; i3 < 3 ; i3 ++) { if (i3 != i2 && i3 != i1) { mT_leptons.push_back(chosen_leptons[i3]); if ( fabs(91.0 - OSSF_inv_mass) < 20. ) mTW_leptons.push_back(chosen_leptons[i3]); } } } else { mT_leptons.push_back(chosen_leptons[0]); mTW_leptons.push_back(chosen_leptons[0]); } } } sortByPt(mT_leptons); sortByPt(mTW_leptons); double mT = sqrt(2*pTmiss.pT()/GeV*mT_leptons[0].pT()/GeV*(1-cos(pTmiss.phi()-mT_leptons[0].phi()))); double mTW = sqrt(2*pTmiss.pT()/GeV*mTW_leptons[0].pT()/GeV*(1-cos(pTmiss.phi()-mTW_leptons[0].phi()))); // Calculate Min pT variable double min_pT = chosen_leptons[2].pT()/GeV; // 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); + _h_e_n->fill(recon_e.size()); + _h_mu_n->fill(recon_mu.size()); + _h_tau_n->fill(recon_tau.size()); // Calculate HTjets variable double HTjets = 0.; for (const Jet& jet : recon_jets) HTjets += jet.pT()/GeV; // Calculate meff variable double meff = eTmiss + HTjets; Particles all_leptons; for (const Particle& e : recon_e ) { meff += e.pT()/GeV; all_leptons.push_back( e ); } for (const Particle& mu : recon_mu) { meff += mu.pT()/GeV; all_leptons.push_back( mu ); } for (const Particle& tau : recon_tau) { meff += tau.pT()/GeV; all_leptons.push_back( tau ); } // Fill histograms 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); - _h_min_pT_all->fill(min_pT,weight); - _h_mT_all->fill(mT,weight); + _h_HTlep_all->fill(HTlep); + _h_HTjets_all->fill(HTjets); + _h_MET_all->fill(eTmiss); + _h_Meff_all->fill(meff); + _h_min_pT_all->fill(min_pT); + _h_mT_all->fill(mT); // Determine signal region (3l / 2ltau , onZ / offZ OSSF / offZ no-OSSF) // 3l vs. 2ltau 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) { bool OSSF = isOSSF(chosen_leptons); if (OSSF) basic_signal_region += "offZ_OSSF"; else basic_signal_region += "offZ_noOSSF"; } // Check in which signal regions this event falls and adjust event counters // INFO: The b-jet signal regions of the paper are not included in this Rivet implementation - fillEventCountsPerSR(basic_signal_region,onZ,HTlep,eTmiss,HTjets,meff,min_pT,mTW,weight); + fillEventCountsPerSR(basic_signal_region,onZ,HTlep,eTmiss,HTjets,meff,min_pT,mTW); } /// Normalise histograms etc., after the run void finalize() { // Normalize to an integrated luminosity of 1 fb-1 double norm = crossSection()/femtobarn/sumOfWeights(); 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 (size_t i = 0; i < _signal_regions.size(); i++) { - double signal_events = _eventCountsPerSR[_signal_regions[i]] * norm; + double signal_events = _eventCountsPerSR[_signal_regions[i]]->val() * norm; // Use expected upper limits to find best signal region: double UL95 = getUpperLimit(_signal_regions[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 signal_events_best_SR = _eventCountsPerSR[best_signal_region]->val() * norm; double exp_UL_best_SR = getUpperLimit(best_signal_region, false); double obs_UL_best_SR = getUpperLimit(best_signal_region, true); // Print out result cout << "----------------------------------------------------------------------------------------" << endl; cout << "Number of total events: " << sumOfWeights() << 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 << "Efficiency*Acceptance: " << _eventCountsPerSR[best_signal_region]->val()/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% C.L." << 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% C.L." << 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 << "INFO: The b-jet signal regions of the paper are not included in this Rivet implementation." << endl; 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_min_pT_all, norm); scale(_h_mT_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, norm); } } /// Helper functions //@{ /// Function giving a list of all signal regions vector getSignalRegions() { // List of basic signal regions vector basic_signal_regions; basic_signal_regions.push_back("3l_offZ_OSSF"); basic_signal_regions.push_back("3l_offZ_noOSSF"); basic_signal_regions.push_back("3l_onZ"); basic_signal_regions.push_back("2ltau_offZ_OSSF"); basic_signal_regions.push_back("2ltau_offZ_noOSSF"); basic_signal_regions.push_back("2ltau_onZ"); // List of kinematic variables 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"); kinematic_variables.push_back("MeffMt"); kinematic_variables.push_back("MinPt"); vector signal_regions; // Loop over all kinematic variables and basic signal regions for (size_t i0 = 0; i0 < kinematic_variables.size(); i0++) { for (size_t i1 = 0; i1 < basic_signal_regions.size(); i1++) { // Is signal region onZ? int onZ = (basic_signal_regions[i1].find("onZ") != string::npos) ? 1 : 0; // Get cut values for this kinematic variable vector cut_values = getCutsPerSignalRegion(kinematic_variables[i0], onZ); // Loop over all cut values for (size_t i2 = 0; i2 < cut_values.size(); i2++) { // Push signal region into vector signal_regions.push_back( kinematic_variables[i0] + "_" + basic_signal_regions[i1] + "_cut_" + toString(cut_values[i2]) ); } } } return signal_regions; } /// Function giving all cut values per kinematic variable vector getCutsPerSignalRegion(const string& signal_region, int onZ = 0) { vector cutValues; // Cut values for HTlep if (signal_region.compare("HTlep") == 0) { cutValues.push_back(0); cutValues.push_back(200); cutValues.push_back(500); cutValues.push_back(800); } // Cut values for MinPt else if (signal_region.compare("MinPt") == 0) { cutValues.push_back(0); cutValues.push_back(50); cutValues.push_back(100); cutValues.push_back(150); } // Cut values for METStrong (HTjets > 150 GeV) and METWeak (HTjets < 150 GeV) else if (signal_region.compare("METStrong") == 0 || signal_region.compare("METWeak") == 0) { cutValues.push_back(0); cutValues.push_back(100); cutValues.push_back(200); cutValues.push_back(300); } // Cut values for Meff if (signal_region.compare("Meff") == 0) { cutValues.push_back(0); cutValues.push_back(600); cutValues.push_back(1000); cutValues.push_back(1500); } // Cut values for MeffStrong (MET > 100 GeV) if ((signal_region.compare("MeffStrong") == 0 || signal_region.compare("MeffMt") == 0) && onZ ==1) { cutValues.push_back(0); cutValues.push_back(600); cutValues.push_back(1200); } return cutValues; } /// function fills map _eventCountsPerSR by looping over all signal regions /// and looking if the event falls into this signal region void fillEventCountsPerSR(const string& basic_signal_region, int onZ, double HTlep, double eTmiss, double HTjets, - double meff, double min_pT, double mTW, - double weight) { + double meff, double min_pT, double mTW) { // Get cut values for HTlep, loop over them and add event if cut is passed vector cut_values = getCutsPerSignalRegion("HTlep", onZ); for (size_t i = 0; i < cut_values.size(); i++) { if (HTlep > cut_values[i]) - _eventCountsPerSR[("HTlep_" + basic_signal_region + "_cut_" + toString(cut_values[i]))] += weight; + _eventCountsPerSR[("HTlep_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill(); } // Get cut values for MinPt, loop over them and add event if cut is passed cut_values = getCutsPerSignalRegion("MinPt", onZ); for (size_t i = 0; i < cut_values.size(); i++) { if (min_pT > cut_values[i]) - _eventCountsPerSR[("MinPt_" + basic_signal_region + "_cut_" + toString(cut_values[i]))] += weight; + _eventCountsPerSR[("MinPt_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill(); } // Get cut values for METStrong, loop over them and add event if cut is passed cut_values = getCutsPerSignalRegion("METStrong", onZ); for (size_t i = 0; i < cut_values.size(); i++) { if (eTmiss > cut_values[i] && HTjets > 150.) - _eventCountsPerSR[("METStrong_" + basic_signal_region + "_cut_" + toString(cut_values[i]))] += weight; + _eventCountsPerSR[("METStrong_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill(); } // Get cut values for METWeak, loop over them and add event if cut is passed cut_values = getCutsPerSignalRegion("METWeak", onZ); for (size_t i = 0; i < cut_values.size(); i++) { if (eTmiss > cut_values[i] && HTjets <= 150.) - _eventCountsPerSR[("METWeak_" + basic_signal_region + "_cut_" + toString(cut_values[i]))] += weight; + _eventCountsPerSR[("METWeak_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill(); } // Get cut values for Meff, loop over them and add event if cut is passed cut_values = getCutsPerSignalRegion("Meff", onZ); for (size_t i = 0; i < cut_values.size(); i++) { if (meff > cut_values[i]) - _eventCountsPerSR[("Meff_" + basic_signal_region + "_cut_" + toString(cut_values[i]))] += weight; + _eventCountsPerSR[("Meff_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill(); } // Get cut values for MeffStrong, loop over them and add event if cut is passed cut_values = getCutsPerSignalRegion("MeffStrong", onZ); for (size_t i = 0; i < cut_values.size(); i++) { if (meff > cut_values[i] && eTmiss > 100.) - _eventCountsPerSR[("MeffStrong_" + basic_signal_region + "_cut_" + toString(cut_values[i]))] += weight; + _eventCountsPerSR[("MeffStrong_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill(); } // Get cut values for MeffMt, loop over them and add event if cut is passed cut_values = getCutsPerSignalRegion("MeffMt", onZ); for (size_t i = 0; i < cut_values.size(); i++) { if (meff > cut_values[i] && mTW > 100. && onZ == 1) - _eventCountsPerSR[("MeffMt_" + basic_signal_region + "_cut_" + toString(cut_values[i]))] += weight; + _eventCountsPerSR[("MeffMt_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill(); } } /// Function returning 4-momentum of daughter-particle if it is a tau neutrino FourMomentum get_tau_neutrino_momentum(const Particle& p) { assert(p.abspid() == PID::TAU); const GenVertex* dv = p.genParticle()->end_vertex(); assert(dv != NULL); // Loop over all daughter particles for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin(); pp != dv->particles_out_const_end(); ++pp) { if (abs((*pp)->pdg_id()) == PID::NU_TAU) return FourMomentum((*pp)->momentum()); } return FourMomentum(); } /// Function calculating the prong number of taus void get_prong_number(const GenParticle* p, unsigned int& nprong, bool& lep_decaying_tau) { assert(p != NULL); const GenVertex* dv = p->end_vertex(); assert(dv != NULL); for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin(); pp != dv->particles_out_const_end(); ++pp) { // If they have status 1 and are charged they will produce a track and the prong number is +1 if ((*pp)->status() == 1 ) { const int id = (*pp)->pdg_id(); if (Rivet::PID::charge(id) != 0 ) ++nprong; // Check if tau decays leptonically if (( abs(id) == PID::ELECTRON || abs(id) == PID::MUON || abs(id) == PID::TAU ) && abs(p->pdg_id()) == PID::TAU) 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); } } } /// Function giving fiducial lepton efficiency double apply_reco_eff(int flavor, const Particle& p) { double pt = p.pT()/GeV; double eta = p.eta(); double eff = 0.; if (flavor == 11) { // weight prompt electron -- now including data/MC ID SF in eff. double avgrate = 0.685; const static double wz_ele[] = {0.0256,0.522,0.607,0.654,0.708,0.737,0.761,0.784,0.815,0.835,0.851,0.841,0.898}; // double ewz_ele[] = {0.000257,0.00492,0.00524,0.00519,0.00396,0.00449,0.00538,0.00513,0.00773,0.00753,0.0209,0.0964,0.259}; int ibin = 0; if(pt > 10 && pt < 15) ibin = 0; if(pt > 15 && pt < 20) ibin = 1; if(pt > 20 && pt < 25) ibin = 2; if(pt > 25 && pt < 30) ibin = 3; if(pt > 30 && pt < 40) ibin = 4; if(pt > 40 && pt < 50) ibin = 5; if(pt > 50 && pt < 60) ibin = 6; if(pt > 60 && pt < 80) ibin = 7; if(pt > 80 && pt < 100) ibin = 8; if(pt > 100 && pt < 200) ibin = 9; if(pt > 200 && pt < 400) ibin = 10; if(pt > 400 && pt < 600) ibin = 11; if(pt > 600) ibin = 12; double eff_pt = 0.; eff_pt = wz_ele[ibin]; eta = fabs(eta); const static double wz_ele_eta[] = {0.65,0.714,0.722,0.689,0.635,0.615}; // double ewz_ele_eta[] = {0.00642,0.00355,0.00335,0.004,0.00368,0.00422}; ibin = 0; if(eta > 0 && eta < 0.1) ibin = 0; if(eta > 0.1 && eta < 0.5) ibin = 1; if(eta > 0.5 && eta < 1.0) ibin = 2; if(eta > 1.0 && eta < 1.5) ibin = 3; if(eta > 1.5 && eta < 2.0) ibin = 4; if(eta > 2.0 && eta < 2.5) ibin = 5; double eff_eta = 0.; eff_eta = wz_ele_eta[ibin]; eff = (eff_pt * eff_eta) / avgrate; } if (flavor == 12) { // weight electron from tau double avgrate = 0.476; const static double wz_ele[] = {0.00855,0.409,0.442,0.55,0.632,0.616,0.615,0.642,0.72,0.617}; // double ewz_ele[] = {0.000573,0.0291,0.0366,0.0352,0.0363,0.0474,0.0628,0.0709,0.125,0.109}; int ibin = 0; if(pt > 10 && pt < 15) ibin = 0; if(pt > 15 && pt < 20) ibin = 1; if(pt > 20 && pt < 25) ibin = 2; if(pt > 25 && pt < 30) ibin = 3; if(pt > 30 && pt < 40) ibin = 4; if(pt > 40 && pt < 50) ibin = 5; if(pt > 50 && pt < 60) ibin = 6; if(pt > 60 && pt < 80) ibin = 7; if(pt > 80 && pt < 100) ibin = 8; if(pt > 100) ibin = 9; double eff_pt = 0.; eff_pt = wz_ele[ibin]; eta = fabs(eta); const static double wz_ele_eta[] = {0.546,0.5,0.513,0.421,0.47,0.433}; //double ewz_ele_eta[] = {0.0566,0.0257,0.0263,0.0263,0.0303,0.0321}; ibin = 0; if(eta > 0 && eta < 0.1) ibin = 0; if(eta > 0.1 && eta < 0.5) ibin = 1; if(eta > 0.5 && eta < 1.0) ibin = 2; if(eta > 1.0 && eta < 1.5) ibin = 3; if(eta > 1.5 && eta < 2.0) ibin = 4; if(eta > 2.0 && eta < 2.5) ibin = 5; double eff_eta = 0.; eff_eta = wz_ele_eta[ibin]; eff = (eff_pt * eff_eta) / avgrate; } if (flavor == 13) { // weight prompt muon int ibin = 0; if(pt > 10 && pt < 15) ibin = 0; if(pt > 15 && pt < 20) ibin = 1; if(pt > 20 && pt < 25) ibin = 2; if(pt > 25 && pt < 30) ibin = 3; if(pt > 30 && pt < 40) ibin = 4; if(pt > 40 && pt < 50) ibin = 5; if(pt > 50 && pt < 60) ibin = 6; if(pt > 60 && pt < 80) ibin = 7; if(pt > 80 && pt < 100) ibin = 8; if(pt > 100 && pt < 200) ibin = 9; if(pt > 200 && pt < 400) ibin = 10; if(pt > 400) ibin = 11; if(fabs(eta) < 0.1) { const static double wz_mu[] = {0.00705,0.402,0.478,0.49,0.492,0.499,0.527,0.512,0.53,0.528,0.465,0.465}; //double ewz_mu[] = {0.000298,0.0154,0.017,0.0158,0.0114,0.0123,0.0155,0.0133,0.0196,0.0182,0.0414,0.0414}; double eff_pt = 0.; eff_pt = wz_mu[ibin]; eff = eff_pt; } if(fabs(eta) > 0.1) { const static double wz_mu[] = {0.0224,0.839,0.887,0.91,0.919,0.923,0.925,0.925,0.922,0.918,0.884,0.834}; //double ewz_mu[] = {0.000213,0.00753,0.0074,0.007,0.00496,0.00534,0.00632,0.00583,0.00849,0.00804,0.0224,0.0963}; double eff_pt = 0.; eff_pt = wz_mu[ibin]; eff = eff_pt; } } if (flavor == 14) { // weight muon from tau int ibin = 0; if(pt > 10 && pt < 15) ibin = 0; if(pt > 15 && pt < 20) ibin = 1; if(pt > 20 && pt < 25) ibin = 2; if(pt > 25 && pt < 30) ibin = 3; if(pt > 30 && pt < 40) ibin = 4; if(pt > 40 && pt < 50) ibin = 5; if(pt > 50 && pt < 60) ibin = 6; if(pt > 60 && pt < 80) ibin = 7; if(pt > 80 && pt < 100) ibin = 8; if(pt > 100) ibin = 9; if(fabs(eta) < 0.1) { const static double wz_mu[] = {0.0,0.664,0.124,0.133,0.527,0.283,0.495,0.25,0.5,0.331}; //double ewz_mu[] = {0.0,0.192,0.0437,0.0343,0.128,0.107,0.202,0.125,0.25,0.191}; double eff_pt = 0.; eff_pt = wz_mu[ibin]; eff = eff_pt; } if(fabs(eta) > 0.1) { const static double wz_mu[] = {0.0,0.617,0.655,0.676,0.705,0.738,0.712,0.783,0.646,0.745}; //double ewz_mu[] = {0.0,0.043,0.0564,0.0448,0.0405,0.0576,0.065,0.0825,0.102,0.132}; double eff_pt = 0.; eff_pt = wz_mu[ibin]; eff = eff_pt; } } if (flavor == 15) { // weight hadronic tau 1p double avgrate = 0.16; const static double wz_tau1p[] = {0.0,0.0311,0.148,0.229,0.217,0.292,0.245,0.307,0.227,0.277}; //double ewz_tau1p[] = {0.0,0.00211,0.0117,0.0179,0.0134,0.0248,0.0264,0.0322,0.0331,0.0427}; int ibin = 0; if(pt > 10 && pt < 15) ibin = 0; if(pt > 15 && pt < 20) ibin = 1; if(pt > 20 && pt < 25) ibin = 2; if(pt > 25 && pt < 30) ibin = 3; if(pt > 30 && pt < 40) ibin = 4; if(pt > 40 && pt < 50) ibin = 5; if(pt > 50 && pt < 60) ibin = 6; if(pt > 60 && pt < 80) ibin = 7; if(pt > 80 && pt < 100) ibin = 8; if(pt > 100) ibin = 9; double eff_pt = 0.; eff_pt = wz_tau1p[ibin]; const static double wz_tau1p_eta[] = {0.166,0.15,0.188,0.175,0.142,0.109}; //double ewz_tau1p_eta[] ={0.0166,0.00853,0.0097,0.00985,0.00949,0.00842}; ibin = 0; if(eta > 0.0 && eta < 0.1) ibin = 0; if(eta > 0.1 && eta < 0.5) ibin = 1; if(eta > 0.5 && eta < 1.0) ibin = 2; if(eta > 1.0 && eta < 1.5) ibin = 3; if(eta > 1.5 && eta < 2.0) ibin = 4; if(eta > 2.0 && eta < 2.5) ibin = 5; double eff_eta = 0.; eff_eta = wz_tau1p_eta[ibin]; eff = (eff_pt * eff_eta) / avgrate; } return eff; } /// Function giving observed and expected upper limits (on the visible cross-section) double getUpperLimit(const string& signal_region, bool observed) { map upperLimitsObserved; map upperLimitsExpected; upperLimitsObserved["HTlep_3l_offZ_OSSF_cut_0"] = 2.435; upperLimitsObserved["HTlep_3l_offZ_OSSF_cut_200"] = 0.704; upperLimitsObserved["HTlep_3l_offZ_OSSF_cut_500"] = 0.182; upperLimitsObserved["HTlep_3l_offZ_OSSF_cut_800"] = 0.147; upperLimitsObserved["HTlep_2ltau_offZ_OSSF_cut_0"] = 13.901; upperLimitsObserved["HTlep_2ltau_offZ_OSSF_cut_200"] = 1.677; upperLimitsObserved["HTlep_2ltau_offZ_OSSF_cut_500"] = 0.141; upperLimitsObserved["HTlep_2ltau_offZ_OSSF_cut_800"] = 0.155; upperLimitsObserved["HTlep_3l_offZ_noOSSF_cut_0"] = 1.054; upperLimitsObserved["HTlep_3l_offZ_noOSSF_cut_200"] = 0.341; upperLimitsObserved["HTlep_3l_offZ_noOSSF_cut_500"] = 0.221; upperLimitsObserved["HTlep_3l_offZ_noOSSF_cut_800"] = 0.140; upperLimitsObserved["HTlep_2ltau_offZ_noOSSF_cut_0"] = 4.276; upperLimitsObserved["HTlep_2ltau_offZ_noOSSF_cut_200"] = 0.413; upperLimitsObserved["HTlep_2ltau_offZ_noOSSF_cut_500"] = 0.138; upperLimitsObserved["HTlep_2ltau_offZ_noOSSF_cut_800"] = 0.150; upperLimitsObserved["HTlep_3l_onZ_cut_0"] = 29.804; upperLimitsObserved["HTlep_3l_onZ_cut_200"] = 3.579; upperLimitsObserved["HTlep_3l_onZ_cut_500"] = 0.466; upperLimitsObserved["HTlep_3l_onZ_cut_800"] = 0.298; upperLimitsObserved["HTlep_2ltau_onZ_cut_0"] = 205.091; upperLimitsObserved["HTlep_2ltau_onZ_cut_200"] = 3.141; upperLimitsObserved["HTlep_2ltau_onZ_cut_500"] = 0.290; upperLimitsObserved["HTlep_2ltau_onZ_cut_800"] = 0.157; upperLimitsObserved["METStrong_3l_offZ_OSSF_cut_0"] = 1.111; upperLimitsObserved["METStrong_3l_offZ_OSSF_cut_100"] = 0.354; upperLimitsObserved["METStrong_3l_offZ_OSSF_cut_200"] = 0.236; upperLimitsObserved["METStrong_3l_offZ_OSSF_cut_300"] = 0.150; upperLimitsObserved["METStrong_2ltau_offZ_OSSF_cut_0"] = 1.881; upperLimitsObserved["METStrong_2ltau_offZ_OSSF_cut_100"] = 0.406; upperLimitsObserved["METStrong_2ltau_offZ_OSSF_cut_200"] = 0.194; upperLimitsObserved["METStrong_2ltau_offZ_OSSF_cut_300"] = 0.134; upperLimitsObserved["METStrong_3l_offZ_noOSSF_cut_0"] = 0.770; upperLimitsObserved["METStrong_3l_offZ_noOSSF_cut_100"] = 0.295; upperLimitsObserved["METStrong_3l_offZ_noOSSF_cut_200"] = 0.149; upperLimitsObserved["METStrong_3l_offZ_noOSSF_cut_300"] = 0.140; upperLimitsObserved["METStrong_2ltau_offZ_noOSSF_cut_0"] = 2.003; upperLimitsObserved["METStrong_2ltau_offZ_noOSSF_cut_100"] = 0.806; upperLimitsObserved["METStrong_2ltau_offZ_noOSSF_cut_200"] = 0.227; upperLimitsObserved["METStrong_2ltau_offZ_noOSSF_cut_300"] = 0.138; upperLimitsObserved["METStrong_3l_onZ_cut_0"] = 6.383; upperLimitsObserved["METStrong_3l_onZ_cut_100"] = 0.959; upperLimitsObserved["METStrong_3l_onZ_cut_200"] = 0.549; upperLimitsObserved["METStrong_3l_onZ_cut_300"] = 0.182; upperLimitsObserved["METStrong_2ltau_onZ_cut_0"] = 10.658; upperLimitsObserved["METStrong_2ltau_onZ_cut_100"] = 0.637; upperLimitsObserved["METStrong_2ltau_onZ_cut_200"] = 0.291; upperLimitsObserved["METStrong_2ltau_onZ_cut_300"] = 0.227; upperLimitsObserved["METWeak_3l_offZ_OSSF_cut_0"] = 1.802; upperLimitsObserved["METWeak_3l_offZ_OSSF_cut_100"] = 0.344; upperLimitsObserved["METWeak_3l_offZ_OSSF_cut_200"] = 0.189; upperLimitsObserved["METWeak_3l_offZ_OSSF_cut_300"] = 0.148; upperLimitsObserved["METWeak_2ltau_offZ_OSSF_cut_0"] = 12.321; upperLimitsObserved["METWeak_2ltau_offZ_OSSF_cut_100"] = 0.430; upperLimitsObserved["METWeak_2ltau_offZ_OSSF_cut_200"] = 0.137; upperLimitsObserved["METWeak_2ltau_offZ_OSSF_cut_300"] = 0.134; upperLimitsObserved["METWeak_3l_offZ_noOSSF_cut_0"] = 0.562; upperLimitsObserved["METWeak_3l_offZ_noOSSF_cut_100"] = 0.153; upperLimitsObserved["METWeak_3l_offZ_noOSSF_cut_200"] = 0.154; upperLimitsObserved["METWeak_3l_offZ_noOSSF_cut_300"] = 0.141; upperLimitsObserved["METWeak_2ltau_offZ_noOSSF_cut_0"] = 2.475; upperLimitsObserved["METWeak_2ltau_offZ_noOSSF_cut_100"] = 0.244; upperLimitsObserved["METWeak_2ltau_offZ_noOSSF_cut_200"] = 0.141; upperLimitsObserved["METWeak_2ltau_offZ_noOSSF_cut_300"] = 0.142; upperLimitsObserved["METWeak_3l_onZ_cut_0"] = 24.769; upperLimitsObserved["METWeak_3l_onZ_cut_100"] = 0.690; upperLimitsObserved["METWeak_3l_onZ_cut_200"] = 0.198; upperLimitsObserved["METWeak_3l_onZ_cut_300"] = 0.138; upperLimitsObserved["METWeak_2ltau_onZ_cut_0"] = 194.360; upperLimitsObserved["METWeak_2ltau_onZ_cut_100"] = 0.287; upperLimitsObserved["METWeak_2ltau_onZ_cut_200"] = 0.144; upperLimitsObserved["METWeak_2ltau_onZ_cut_300"] = 0.130; upperLimitsObserved["Meff_3l_offZ_OSSF_cut_0"] = 2.435; upperLimitsObserved["Meff_3l_offZ_OSSF_cut_600"] = 0.487; upperLimitsObserved["Meff_3l_offZ_OSSF_cut_1000"] = 0.156; upperLimitsObserved["Meff_3l_offZ_OSSF_cut_1500"] = 0.140; upperLimitsObserved["Meff_2ltau_offZ_OSSF_cut_0"] = 13.901; upperLimitsObserved["Meff_2ltau_offZ_OSSF_cut_600"] = 0.687; upperLimitsObserved["Meff_2ltau_offZ_OSSF_cut_1000"] = 0.224; upperLimitsObserved["Meff_2ltau_offZ_OSSF_cut_1500"] = 0.155; upperLimitsObserved["Meff_3l_offZ_noOSSF_cut_0"] = 1.054; upperLimitsObserved["Meff_3l_offZ_noOSSF_cut_600"] = 0.249; upperLimitsObserved["Meff_3l_offZ_noOSSF_cut_1000"] = 0.194; upperLimitsObserved["Meff_3l_offZ_noOSSF_cut_1500"] = 0.145; upperLimitsObserved["Meff_2ltau_offZ_noOSSF_cut_0"] = 4.276; upperLimitsObserved["Meff_2ltau_offZ_noOSSF_cut_600"] = 0.772; upperLimitsObserved["Meff_2ltau_offZ_noOSSF_cut_1000"] = 0.218; upperLimitsObserved["Meff_2ltau_offZ_noOSSF_cut_1500"] = 0.204; upperLimitsObserved["Meff_3l_onZ_cut_0"] = 29.804; upperLimitsObserved["Meff_3l_onZ_cut_600"] = 2.933; upperLimitsObserved["Meff_3l_onZ_cut_1000"] = 0.912; upperLimitsObserved["Meff_3l_onZ_cut_1500"] = 0.225; upperLimitsObserved["Meff_2ltau_onZ_cut_0"] = 205.091; upperLimitsObserved["Meff_2ltau_onZ_cut_600"] = 1.486; upperLimitsObserved["Meff_2ltau_onZ_cut_1000"] = 0.641; upperLimitsObserved["Meff_2ltau_onZ_cut_1500"] = 0.204; upperLimitsObserved["MeffStrong_3l_offZ_OSSF_cut_0"] = 0.479; upperLimitsObserved["MeffStrong_3l_offZ_OSSF_cut_600"] = 0.353; upperLimitsObserved["MeffStrong_3l_offZ_OSSF_cut_1200"] = 0.187; upperLimitsObserved["MeffStrong_2ltau_offZ_OSSF_cut_0"] = 0.617; upperLimitsObserved["MeffStrong_2ltau_offZ_OSSF_cut_600"] = 0.320; upperLimitsObserved["MeffStrong_2ltau_offZ_OSSF_cut_1200"] = 0.281; upperLimitsObserved["MeffStrong_3l_offZ_noOSSF_cut_0"] = 0.408; upperLimitsObserved["MeffStrong_3l_offZ_noOSSF_cut_600"] = 0.240; upperLimitsObserved["MeffStrong_3l_offZ_noOSSF_cut_1200"] = 0.150; upperLimitsObserved["MeffStrong_2ltau_offZ_noOSSF_cut_0"] = 0.774; upperLimitsObserved["MeffStrong_2ltau_offZ_noOSSF_cut_600"] = 0.417; upperLimitsObserved["MeffStrong_2ltau_offZ_noOSSF_cut_1200"] = 0.266; upperLimitsObserved["MeffStrong_3l_onZ_cut_0"] = 1.208; upperLimitsObserved["MeffStrong_3l_onZ_cut_600"] = 0.837; upperLimitsObserved["MeffStrong_3l_onZ_cut_1200"] = 0.269; upperLimitsObserved["MeffStrong_2ltau_onZ_cut_0"] = 0.605; upperLimitsObserved["MeffStrong_2ltau_onZ_cut_600"] = 0.420; upperLimitsObserved["MeffStrong_2ltau_onZ_cut_1200"] = 0.141; upperLimitsObserved["MeffMt_3l_onZ_cut_0"] = 1.832; upperLimitsObserved["MeffMt_3l_onZ_cut_600"] = 0.862; upperLimitsObserved["MeffMt_3l_onZ_cut_1200"] = 0.222; upperLimitsObserved["MeffMt_2ltau_onZ_cut_0"] = 1.309; upperLimitsObserved["MeffMt_2ltau_onZ_cut_600"] = 0.481; upperLimitsObserved["MeffMt_2ltau_onZ_cut_1200"] = 0.146; upperLimitsObserved["MinPt_3l_offZ_OSSF_cut_0"] = 2.435; upperLimitsObserved["MinPt_3l_offZ_OSSF_cut_50"] = 0.500; upperLimitsObserved["MinPt_3l_offZ_OSSF_cut_100"] = 0.203; upperLimitsObserved["MinPt_3l_offZ_OSSF_cut_150"] = 0.128; upperLimitsObserved["MinPt_2ltau_offZ_OSSF_cut_0"] = 13.901; upperLimitsObserved["MinPt_2ltau_offZ_OSSF_cut_50"] = 0.859; upperLimitsObserved["MinPt_2ltau_offZ_OSSF_cut_100"] = 0.158; upperLimitsObserved["MinPt_2ltau_offZ_OSSF_cut_150"] = 0.155; upperLimitsObserved["MinPt_3l_offZ_noOSSF_cut_0"] = 1.054; upperLimitsObserved["MinPt_3l_offZ_noOSSF_cut_50"] = 0.295; upperLimitsObserved["MinPt_3l_offZ_noOSSF_cut_100"] = 0.148; upperLimitsObserved["MinPt_3l_offZ_noOSSF_cut_150"] = 0.137; upperLimitsObserved["MinPt_2ltau_offZ_noOSSF_cut_0"] = 4.276; upperLimitsObserved["MinPt_2ltau_offZ_noOSSF_cut_50"] = 0.314; upperLimitsObserved["MinPt_2ltau_offZ_noOSSF_cut_100"] = 0.134; upperLimitsObserved["MinPt_2ltau_offZ_noOSSF_cut_150"] = 0.140; upperLimitsObserved["MinPt_3l_onZ_cut_0"] = 29.804; upperLimitsObserved["MinPt_3l_onZ_cut_50"] = 1.767; upperLimitsObserved["MinPt_3l_onZ_cut_100"] = 0.690; upperLimitsObserved["MinPt_3l_onZ_cut_150"] = 0.301; upperLimitsObserved["MinPt_2ltau_onZ_cut_0"] = 205.091; upperLimitsObserved["MinPt_2ltau_onZ_cut_50"] = 1.050; upperLimitsObserved["MinPt_2ltau_onZ_cut_100"] = 0.155; upperLimitsObserved["MinPt_2ltau_onZ_cut_150"] = 0.146; upperLimitsObserved["nbtag_3l_offZ_OSSF_cut_0"] = 2.435; upperLimitsObserved["nbtag_3l_offZ_OSSF_cut_1"] = 0.865; upperLimitsObserved["nbtag_3l_offZ_OSSF_cut_2"] = 0.474; upperLimitsObserved["nbtag_2ltau_offZ_OSSF_cut_0"] = 13.901; upperLimitsObserved["nbtag_2ltau_offZ_OSSF_cut_1"] = 1.566; upperLimitsObserved["nbtag_2ltau_offZ_OSSF_cut_2"] = 0.426; upperLimitsObserved["nbtag_3l_offZ_noOSSF_cut_0"] = 1.054; upperLimitsObserved["nbtag_3l_offZ_noOSSF_cut_1"] = 0.643; upperLimitsObserved["nbtag_3l_offZ_noOSSF_cut_2"] = 0.321; upperLimitsObserved["nbtag_2ltau_offZ_noOSSF_cut_0"] = 4.276; upperLimitsObserved["nbtag_2ltau_offZ_noOSSF_cut_1"] = 2.435; upperLimitsObserved["nbtag_2ltau_offZ_noOSSF_cut_2"] = 1.073; upperLimitsObserved["nbtag_3l_onZ_cut_0"] = 29.804; upperLimitsObserved["nbtag_3l_onZ_cut_1"] = 3.908; upperLimitsObserved["nbtag_3l_onZ_cut_2"] = 0.704; upperLimitsObserved["nbtag_2ltau_onZ_cut_0"] = 205.091; upperLimitsObserved["nbtag_2ltau_onZ_cut_1"] = 9.377; upperLimitsObserved["nbtag_2ltau_onZ_cut_2"] = 0.657; upperLimitsExpected["HTlep_3l_offZ_OSSF_cut_0"] = 2.893; upperLimitsExpected["HTlep_3l_offZ_OSSF_cut_200"] = 1.175; upperLimitsExpected["HTlep_3l_offZ_OSSF_cut_500"] = 0.265; upperLimitsExpected["HTlep_3l_offZ_OSSF_cut_800"] = 0.155; upperLimitsExpected["HTlep_2ltau_offZ_OSSF_cut_0"] = 14.293; upperLimitsExpected["HTlep_2ltau_offZ_OSSF_cut_200"] = 1.803; upperLimitsExpected["HTlep_2ltau_offZ_OSSF_cut_500"] = 0.159; upperLimitsExpected["HTlep_2ltau_offZ_OSSF_cut_800"] = 0.155; upperLimitsExpected["HTlep_3l_offZ_noOSSF_cut_0"] = 0.836; upperLimitsExpected["HTlep_3l_offZ_noOSSF_cut_200"] = 0.340; upperLimitsExpected["HTlep_3l_offZ_noOSSF_cut_500"] = 0.218; upperLimitsExpected["HTlep_3l_offZ_noOSSF_cut_800"] = 0.140; upperLimitsExpected["HTlep_2ltau_offZ_noOSSF_cut_0"] = 4.132; upperLimitsExpected["HTlep_2ltau_offZ_noOSSF_cut_200"] = 0.599; upperLimitsExpected["HTlep_2ltau_offZ_noOSSF_cut_500"] = 0.146; upperLimitsExpected["HTlep_2ltau_offZ_noOSSF_cut_800"] = 0.148; upperLimitsExpected["HTlep_3l_onZ_cut_0"] = 32.181; upperLimitsExpected["HTlep_3l_onZ_cut_200"] = 4.879; upperLimitsExpected["HTlep_3l_onZ_cut_500"] = 0.473; upperLimitsExpected["HTlep_3l_onZ_cut_800"] = 0.266; upperLimitsExpected["HTlep_2ltau_onZ_cut_0"] = 217.801; upperLimitsExpected["HTlep_2ltau_onZ_cut_200"] = 3.676; upperLimitsExpected["HTlep_2ltau_onZ_cut_500"] = 0.235; upperLimitsExpected["HTlep_2ltau_onZ_cut_800"] = 0.150; upperLimitsExpected["METStrong_3l_offZ_OSSF_cut_0"] = 1.196; upperLimitsExpected["METStrong_3l_offZ_OSSF_cut_100"] = 0.423; upperLimitsExpected["METStrong_3l_offZ_OSSF_cut_200"] = 0.208; upperLimitsExpected["METStrong_3l_offZ_OSSF_cut_300"] = 0.158; upperLimitsExpected["METStrong_2ltau_offZ_OSSF_cut_0"] = 2.158; upperLimitsExpected["METStrong_2ltau_offZ_OSSF_cut_100"] = 0.461; upperLimitsExpected["METStrong_2ltau_offZ_OSSF_cut_200"] = 0.186; upperLimitsExpected["METStrong_2ltau_offZ_OSSF_cut_300"] = 0.138; upperLimitsExpected["METStrong_3l_offZ_noOSSF_cut_0"] = 0.495; upperLimitsExpected["METStrong_3l_offZ_noOSSF_cut_100"] = 0.284; upperLimitsExpected["METStrong_3l_offZ_noOSSF_cut_200"] = 0.150; upperLimitsExpected["METStrong_3l_offZ_noOSSF_cut_300"] = 0.146; upperLimitsExpected["METStrong_2ltau_offZ_noOSSF_cut_0"] = 1.967; upperLimitsExpected["METStrong_2ltau_offZ_noOSSF_cut_100"] = 0.732; upperLimitsExpected["METStrong_2ltau_offZ_noOSSF_cut_200"] = 0.225; upperLimitsExpected["METStrong_2ltau_offZ_noOSSF_cut_300"] = 0.147; upperLimitsExpected["METStrong_3l_onZ_cut_0"] = 7.157; upperLimitsExpected["METStrong_3l_onZ_cut_100"] = 1.342; upperLimitsExpected["METStrong_3l_onZ_cut_200"] = 0.508; upperLimitsExpected["METStrong_3l_onZ_cut_300"] = 0.228; upperLimitsExpected["METStrong_2ltau_onZ_cut_0"] = 12.441; upperLimitsExpected["METStrong_2ltau_onZ_cut_100"] = 0.534; upperLimitsExpected["METStrong_2ltau_onZ_cut_200"] = 0.243; upperLimitsExpected["METStrong_2ltau_onZ_cut_300"] = 0.218; upperLimitsExpected["METWeak_3l_offZ_OSSF_cut_0"] = 2.199; upperLimitsExpected["METWeak_3l_offZ_OSSF_cut_100"] = 0.391; upperLimitsExpected["METWeak_3l_offZ_OSSF_cut_200"] = 0.177; upperLimitsExpected["METWeak_3l_offZ_OSSF_cut_300"] = 0.144; upperLimitsExpected["METWeak_2ltau_offZ_OSSF_cut_0"] = 12.431; upperLimitsExpected["METWeak_2ltau_offZ_OSSF_cut_100"] = 0.358; upperLimitsExpected["METWeak_2ltau_offZ_OSSF_cut_200"] = 0.150; upperLimitsExpected["METWeak_2ltau_offZ_OSSF_cut_300"] = 0.135; upperLimitsExpected["METWeak_3l_offZ_noOSSF_cut_0"] = 0.577; upperLimitsExpected["METWeak_3l_offZ_noOSSF_cut_100"] = 0.214; upperLimitsExpected["METWeak_3l_offZ_noOSSF_cut_200"] = 0.155; upperLimitsExpected["METWeak_3l_offZ_noOSSF_cut_300"] = 0.140; upperLimitsExpected["METWeak_2ltau_offZ_noOSSF_cut_0"] = 2.474; upperLimitsExpected["METWeak_2ltau_offZ_noOSSF_cut_100"] = 0.382; upperLimitsExpected["METWeak_2ltau_offZ_noOSSF_cut_200"] = 0.144; upperLimitsExpected["METWeak_2ltau_offZ_noOSSF_cut_300"] = 0.146; upperLimitsExpected["METWeak_3l_onZ_cut_0"] = 26.305; upperLimitsExpected["METWeak_3l_onZ_cut_100"] = 1.227; upperLimitsExpected["METWeak_3l_onZ_cut_200"] = 0.311; upperLimitsExpected["METWeak_3l_onZ_cut_300"] = 0.188; upperLimitsExpected["METWeak_2ltau_onZ_cut_0"] = 205.198; upperLimitsExpected["METWeak_2ltau_onZ_cut_100"] = 0.399; upperLimitsExpected["METWeak_2ltau_onZ_cut_200"] = 0.166; upperLimitsExpected["METWeak_2ltau_onZ_cut_300"] = 0.140; upperLimitsExpected["Meff_3l_offZ_OSSF_cut_0"] = 2.893; upperLimitsExpected["Meff_3l_offZ_OSSF_cut_600"] = 0.649; upperLimitsExpected["Meff_3l_offZ_OSSF_cut_1000"] = 0.252; upperLimitsExpected["Meff_3l_offZ_OSSF_cut_1500"] = 0.150; upperLimitsExpected["Meff_2ltau_offZ_OSSF_cut_0"] = 14.293; upperLimitsExpected["Meff_2ltau_offZ_OSSF_cut_600"] = 0.657; upperLimitsExpected["Meff_2ltau_offZ_OSSF_cut_1000"] = 0.226; upperLimitsExpected["Meff_2ltau_offZ_OSSF_cut_1500"] = 0.154; upperLimitsExpected["Meff_3l_offZ_noOSSF_cut_0"] = 0.836; upperLimitsExpected["Meff_3l_offZ_noOSSF_cut_600"] = 0.265; upperLimitsExpected["Meff_3l_offZ_noOSSF_cut_1000"] = 0.176; upperLimitsExpected["Meff_3l_offZ_noOSSF_cut_1500"] = 0.146; upperLimitsExpected["Meff_2ltau_offZ_noOSSF_cut_0"] = 4.132; upperLimitsExpected["Meff_2ltau_offZ_noOSSF_cut_600"] = 0.678; upperLimitsExpected["Meff_2ltau_offZ_noOSSF_cut_1000"] = 0.243; upperLimitsExpected["Meff_2ltau_offZ_noOSSF_cut_1500"] = 0.184; upperLimitsExpected["Meff_3l_onZ_cut_0"] = 32.181; upperLimitsExpected["Meff_3l_onZ_cut_600"] = 3.219; upperLimitsExpected["Meff_3l_onZ_cut_1000"] = 0.905; upperLimitsExpected["Meff_3l_onZ_cut_1500"] = 0.261; upperLimitsExpected["Meff_2ltau_onZ_cut_0"] = 217.801; upperLimitsExpected["Meff_2ltau_onZ_cut_600"] = 1.680; upperLimitsExpected["Meff_2ltau_onZ_cut_1000"] = 0.375; upperLimitsExpected["Meff_2ltau_onZ_cut_1500"] = 0.178; upperLimitsExpected["MeffStrong_3l_offZ_OSSF_cut_0"] = 0.571; upperLimitsExpected["MeffStrong_3l_offZ_OSSF_cut_600"] = 0.386; upperLimitsExpected["MeffStrong_3l_offZ_OSSF_cut_1200"] = 0.177; upperLimitsExpected["MeffStrong_2ltau_offZ_OSSF_cut_0"] = 0.605; upperLimitsExpected["MeffStrong_2ltau_offZ_OSSF_cut_600"] = 0.335; upperLimitsExpected["MeffStrong_2ltau_offZ_OSSF_cut_1200"] = 0.249; upperLimitsExpected["MeffStrong_3l_offZ_noOSSF_cut_0"] = 0.373; upperLimitsExpected["MeffStrong_3l_offZ_noOSSF_cut_600"] = 0.223; upperLimitsExpected["MeffStrong_3l_offZ_noOSSF_cut_1200"] = 0.150; upperLimitsExpected["MeffStrong_2ltau_offZ_noOSSF_cut_0"] = 0.873; upperLimitsExpected["MeffStrong_2ltau_offZ_noOSSF_cut_600"] = 0.428; upperLimitsExpected["MeffStrong_2ltau_offZ_noOSSF_cut_1200"] = 0.210; upperLimitsExpected["MeffStrong_3l_onZ_cut_0"] = 2.034; upperLimitsExpected["MeffStrong_3l_onZ_cut_600"] = 1.093; upperLimitsExpected["MeffStrong_3l_onZ_cut_1200"] = 0.293; upperLimitsExpected["MeffStrong_2ltau_onZ_cut_0"] = 0.690; upperLimitsExpected["MeffStrong_2ltau_onZ_cut_600"] = 0.392; upperLimitsExpected["MeffStrong_2ltau_onZ_cut_1200"] = 0.156; upperLimitsExpected["MeffMt_3l_onZ_cut_0"] = 2.483; upperLimitsExpected["MeffMt_3l_onZ_cut_600"] = 0.845; upperLimitsExpected["MeffMt_3l_onZ_cut_1200"] = 0.255; upperLimitsExpected["MeffMt_2ltau_onZ_cut_0"] = 1.448; upperLimitsExpected["MeffMt_2ltau_onZ_cut_600"] = 0.391; upperLimitsExpected["MeffMt_2ltau_onZ_cut_1200"] = 0.146; upperLimitsExpected["MinPt_3l_offZ_OSSF_cut_0"] = 2.893; upperLimitsExpected["MinPt_3l_offZ_OSSF_cut_50"] = 0.703; upperLimitsExpected["MinPt_3l_offZ_OSSF_cut_100"] = 0.207; upperLimitsExpected["MinPt_3l_offZ_OSSF_cut_150"] = 0.143; upperLimitsExpected["MinPt_2ltau_offZ_OSSF_cut_0"] = 14.293; upperLimitsExpected["MinPt_2ltau_offZ_OSSF_cut_50"] = 0.705; upperLimitsExpected["MinPt_2ltau_offZ_OSSF_cut_100"] = 0.149; upperLimitsExpected["MinPt_2ltau_offZ_OSSF_cut_150"] = 0.155; upperLimitsExpected["MinPt_3l_offZ_noOSSF_cut_0"] = 0.836; upperLimitsExpected["MinPt_3l_offZ_noOSSF_cut_50"] = 0.249; upperLimitsExpected["MinPt_3l_offZ_noOSSF_cut_100"] = 0.135; upperLimitsExpected["MinPt_3l_offZ_noOSSF_cut_150"] = 0.136; upperLimitsExpected["MinPt_2ltau_offZ_noOSSF_cut_0"] = 4.132; upperLimitsExpected["MinPt_2ltau_offZ_noOSSF_cut_50"] = 0.339; upperLimitsExpected["MinPt_2ltau_offZ_noOSSF_cut_100"] = 0.149; upperLimitsExpected["MinPt_2ltau_offZ_noOSSF_cut_150"] = 0.145; upperLimitsExpected["MinPt_3l_onZ_cut_0"] = 32.181; upperLimitsExpected["MinPt_3l_onZ_cut_50"] = 2.260; upperLimitsExpected["MinPt_3l_onZ_cut_100"] = 0.438; upperLimitsExpected["MinPt_3l_onZ_cut_150"] = 0.305; upperLimitsExpected["MinPt_2ltau_onZ_cut_0"] = 217.801; upperLimitsExpected["MinPt_2ltau_onZ_cut_50"] = 1.335; upperLimitsExpected["MinPt_2ltau_onZ_cut_100"] = 0.162; upperLimitsExpected["MinPt_2ltau_onZ_cut_150"] = 0.149; upperLimitsExpected["nbtag_3l_offZ_OSSF_cut_0"] = 2.893; upperLimitsExpected["nbtag_3l_offZ_OSSF_cut_1"] = 0.923; upperLimitsExpected["nbtag_3l_offZ_OSSF_cut_2"] = 0.452; upperLimitsExpected["nbtag_2ltau_offZ_OSSF_cut_0"] = 14.293; upperLimitsExpected["nbtag_2ltau_offZ_OSSF_cut_1"] = 1.774; upperLimitsExpected["nbtag_2ltau_offZ_OSSF_cut_2"] = 0.549; upperLimitsExpected["nbtag_3l_offZ_noOSSF_cut_0"] = 0.836; upperLimitsExpected["nbtag_3l_offZ_noOSSF_cut_1"] = 0.594; upperLimitsExpected["nbtag_3l_offZ_noOSSF_cut_2"] = 0.298; upperLimitsExpected["nbtag_2ltau_offZ_noOSSF_cut_0"] = 4.132; upperLimitsExpected["nbtag_2ltau_offZ_noOSSF_cut_1"] = 2.358; upperLimitsExpected["nbtag_2ltau_offZ_noOSSF_cut_2"] = 0.958; upperLimitsExpected["nbtag_3l_onZ_cut_0"] = 32.181; upperLimitsExpected["nbtag_3l_onZ_cut_1"] = 3.868; upperLimitsExpected["nbtag_3l_onZ_cut_2"] = 0.887; upperLimitsExpected["nbtag_2ltau_onZ_cut_0"] = 217.801; upperLimitsExpected["nbtag_2ltau_onZ_cut_1"] = 9.397; upperLimitsExpected["nbtag_2ltau_onZ_cut_2"] = 0.787; if (observed) return upperLimitsObserved[signal_region]; else return upperLimitsExpected[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 Particles& 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 for (const Particle& p1 : particles) { for (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.pid()*p2.pid() == -121 || p1.pid()*p2.pid() == -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) for (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; } /// function checking if two leptons are an OSSF lepton pair and giving out the invariant mass (0 if no OSSF pair) double isOSSF_mass (const Particle& p1, const Particle& p2) { double inv_mass = 0.; // Is particle combination OSSF pair? if ((p1.pid()*p2.pid() == -121 || p1.pid()*p2.pid() == -169)) { // Get invariant mass inv_mass = (p1.momentum() + p2.momentum()).mass()/GeV; } return inv_mass; } /// Function checking if there is an OSSF lepton pair bool isOSSF (const Particles& particles) { for (size_t i1=0 ; i1 < 3 ; i1 ++) { for (size_t i2 = i1+1 ; i2 < 3 ; i2 ++) { if ((particles[i1].pid()*particles[i2].pid() == -121 || particles[i1].pid()*particles[i2].pid() == -169)) { return true; } } } return false; } //@} private: /// Histograms //@{ Histo1DPtr _h_HTlep_all, _h_HTjets_all, _h_MET_all, _h_Meff_all, _h_min_pT_all, _h_mT_all; Histo1DPtr _h_pt_1_3l, _h_pt_2_3l, _h_pt_3_3l, _h_pt_1_2ltau, _h_pt_2_2ltau, _h_pt_3_2ltau; Histo1DPtr _h_e_n, _h_mu_n, _h_tau_n; Histo1DPtr _h_excluded; //@} /// Fiducial efficiencies to model the effects of the ATLAS detector bool _use_fiducial_lepton_efficiency; /// List of signal regions and event counts per signal region vector _signal_regions; - map _eventCountsPerSR; + map _eventCountsPerSR; }; DECLARE_RIVET_PLUGIN(ATLAS_2014_I1327229); } diff --git a/analyses/pluginATLAS/ATLAS_2015_CONF_2015_041.cc b/analyses/pluginATLAS/ATLAS_2015_CONF_2015_041.cc --- a/analyses/pluginATLAS/ATLAS_2015_CONF_2015_041.cc +++ b/analyses/pluginATLAS/ATLAS_2015_CONF_2015_041.cc @@ -1,160 +1,160 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/ZFinder.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/VetoedFinalState.hh" namespace Rivet { - /// Z + jets in pp at 13 TeV + /// Z + jets in pp at 13 TeV /// @note This base class contains a "mode" variable for combined, e, and mu channel derived classes class ATLAS_2015_CONF_2015_041 : public Analysis { public: /// @name Constructors etc. //@{ /// Constructor ATLAS_2015_CONF_2015_041(string name="ATLAS_2015_CONF_2015_041") - : Analysis(name), - _weights(5, 0.0) + : Analysis(name), _weights(5) { // This class uses the combined e+mu mode _mode = 0; } //@} /// Book histograms and initialise projections before the run void init() { const FinalState fs; Cut cuts = (Cuts::pT > 25*GeV) & (Cuts::abseta < 2.5); ZFinder zfinder(fs, cuts, _mode? PID::MUON : PID::ELECTRON, 66*GeV, 116*GeV); declare(zfinder, "zfinder"); // Define veto FS in order to prevent Z-decay products entering the jet algorithm VetoedFinalState had_fs; had_fs.addVetoOnThisFinalState(zfinder); FastJets jets(had_fs, FastJets::ANTIKT, 0.4); jets.useInvisibles(true); declare(jets, "jets"); // individual channels book(_hNjets ,1, 1, _mode + 1); book(_hNjetsRatio ,2, 1, _mode + 1, true); // combination book(_hNjets_comb ,1, 2, _mode + 1); book(_hNjetsRatio_comb ,2, 2, _mode + 1, true); + + for (size_t i = 0; i < 5; i++) + book(_weights[i], (ostringstream("_weights") << i).str()); } /// Perform the per-event analysis void analyze(const Event& event) { - const double weight = 1.0; - const ZFinder& zfinder = apply(event, "zfinder"); const Particles& leptons = zfinder.constituents(); if (leptons.size() != 2) vetoEvent; Jets jets; foreach (Jet j, apply(event, "jets").jetsByPt(Cuts::pT > 30*GeV && Cuts::absrap < 2.5)) { bool keep = true; foreach(const Particle& l, leptons) keep &= deltaR(j, l) > 0.4; if (keep) jets += j; } size_t njets = jets.size(); for(size_t i = 0; i <= njets; ++i) { - _hNjets->fill(i + 0.5, weight); - _hNjets_comb->fill(i + 0.5, weight); + _hNjets->fill(i + 0.5); + _hNjets_comb->fill(i + 0.5); } for (size_t i = 0; i < 5; ++i) { - if (njets >= i) _weights[i] += weight; + if (njets >= i) _weights[i]->fill(); } } /// @name Ratio calculator util functions //@{ /// Calculate the ratio, being careful about div-by-zero double ratio(double a, double b) { return (b != 0) ? a/b : 0; } /// Calculate the ratio error, being careful about div-by-zero double ratio_err(double a, double b) { return (b != 0) ? sqrt(a/b*(1-a/b)/b) : 0; } //@} void finalize() { for (size_t i = 0; i < 4; ++i) { double n = _hNjets->bin(i + 1).sumW(); double dN = _hNjets->bin(i + 1).sumW2(); double d = _hNjets->bin(i).sumW(); double dD = _hNjets->bin(i).sumW2(); double r = safediv(n, d); double e = sqrt( safediv(r * (1 - r), d) ); if ( _hNjets->effNumEntries() != _hNjets->numEntries() ) { // use F. James's approximation for weighted events: e = sqrt( safediv((1 - 2 * r) * dN + r * r * dD, d * d) ); } _hNjetsRatio->point(i).setY(r, e); _hNjetsRatio_comb->point(i).setY(r, e); } scale(_hNjets, crossSectionPerEvent() ); scale(_hNjets_comb, crossSectionPerEvent() ); } //@} protected: size_t _mode; private: - vector _weights; + vector _weights; Scatter2DPtr _hNjetsRatio, _hNjetsRatio_comb; Histo1DPtr _hNjets, _hNjets_comb; }; class ATLAS_2015_CONF_2015_041_EL : public ATLAS_2015_CONF_2015_041 { public: ATLAS_2015_CONF_2015_041_EL() : ATLAS_2015_CONF_2015_041("ATLAS_2015_CONF_2015_041_EL") { _mode = 0; } }; class ATLAS_2015_CONF_2015_041_MU : public ATLAS_2015_CONF_2015_041 { public: ATLAS_2015_CONF_2015_041_MU() : ATLAS_2015_CONF_2015_041("ATLAS_2015_CONF_2015_041_MU") { _mode = 1; } }; DECLARE_RIVET_PLUGIN(ATLAS_2015_CONF_2015_041); DECLARE_RIVET_PLUGIN(ATLAS_2015_CONF_2015_041_EL); DECLARE_RIVET_PLUGIN(ATLAS_2015_CONF_2015_041_MU); } diff --git a/analyses/pluginATLAS/ATLAS_2015_I1345452.cc b/analyses/pluginATLAS/ATLAS_2015_I1345452.cc --- a/analyses/pluginATLAS/ATLAS_2015_I1345452.cc +++ b/analyses/pluginATLAS/ATLAS_2015_I1345452.cc @@ -1,289 +1,287 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { /// @brief ATLAS 7 TeV pseudo-top analysis /// /// @author K .Finelli /// @author A. Saavedra /// @author L. Lan class ATLAS_2015_I1345452 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2015_I1345452); void init() { // Eta ranges Cut eta_full = (Cuts::abseta < 5.0) & (Cuts::pT >= 1.0*MeV); Cut eta_lep = (Cuts::abseta < 2.5); // All final state particles FinalState fs(eta_full); // Get photons to dress leptons IdentifiedFinalState photons(fs); photons.acceptIdPair(PID::PHOTON); // Projection to find the electrons IdentifiedFinalState el_id(fs); el_id.acceptIdPair(PID::ELECTRON); PromptFinalState electrons(el_id); electrons.acceptTauDecays(true); declare(electrons, "electrons"); DressedLeptons dressedelectrons(photons, electrons, 0.1, eta_lep && Cuts::pT > 25*GeV, true); declare(dressedelectrons, "dressedelectrons"); DressedLeptons ewdressedelectrons(photons, electrons, 0.1, eta_full, true); declare(ewdressedelectrons, "ewdressedelectrons"); DressedLeptons vetodressedelectrons(photons, electrons, 0.1, eta_lep && Cuts::pT > 15*GeV, true); declare(vetodressedelectrons, "vetodressedelectrons"); // Projection to find the muons IdentifiedFinalState mu_id(fs); mu_id.acceptIdPair(PID::MUON); PromptFinalState muons(mu_id); muons.acceptTauDecays(true); declare(muons, "muons"); DressedLeptons dressedmuons(photons, muons, 0.1, eta_lep && Cuts::pT > 25*GeV, true); declare(dressedmuons, "dressedmuons"); DressedLeptons ewdressedmuons(photons, muons, 0.1, eta_full, true); declare(ewdressedmuons, "ewdressedmuons"); DressedLeptons vetodressedmuons(photons, muons, 0.1, eta_lep && Cuts::pT > 15*GeV, true); declare(vetodressedmuons, "vetodressedmuons"); // Projection to find neutrinos and produce MET IdentifiedFinalState nu_id; nu_id.acceptNeutrinos(); PromptFinalState neutrinos(nu_id); neutrinos.acceptTauDecays(true); declare(neutrinos, "neutrinos"); // Jet clustering. VetoedFinalState vfs; vfs.addVetoOnThisFinalState(ewdressedelectrons); vfs.addVetoOnThisFinalState(ewdressedmuons); vfs.addVetoOnThisFinalState(neutrinos); FastJets jets(vfs, FastJets::ANTIKT, 0.4); jets.useInvisibles(); declare(jets, "jets"); //pseudotop leptons and hadrons book(_h["ptpseudotophadron_mu"] , 1, 1, 2); book(_h["ptpseudotophadron_el"] , 2, 1, 2); book(_h["absrappseudotophadron_mu"] , 3, 1, 2); book(_h["absrappseudotophadron_el"] , 4, 1, 2); book(_h["ptpseudotoplepton_mu"] , 5, 1, 2); book(_h["ptpseudotoplepton_el"] , 6, 1, 2); book(_h["absrappseudotoplepton_mu"] , 7, 1, 2); book(_h["absrappseudotoplepton_el"] , 8, 1, 2); book(_h["ptttbar_mu"] , 9, 1, 2); book(_h["ptttbar_el"] ,10, 1, 2); book(_h["absrapttbar_mu"] ,11, 1, 2); book(_h["absrapttbar_el"] ,12, 1, 2); book(_h["ttbarmass_mu"] ,13, 1, 2); book(_h["ttbarmass_el"] ,14, 1, 2); book(_h["ptpseudotophadron"] ,15, 1, 2); book(_h["absrappseudotophadron"] ,16, 1, 2); book(_h["ptpseudotoplepton"] ,17, 1, 2); book(_h["absrappseudotoplepton"] ,18, 1, 2); book(_h["ptttbar"] ,19, 1, 2); book(_h["absrapttbar"] ,20, 1, 2); book(_h["ttbarmass"] ,21, 1, 2); } void analyze(const Event& event) { // Get the selected objects, using the projections. _dressedelectrons = apply( event, "dressedelectrons").dressedLeptons(); _vetodressedelectrons = apply( event, "vetodressedelectrons").dressedLeptons(); _dressedmuons = apply( event, "dressedmuons").dressedLeptons(); _vetodressedmuons = apply( event, "vetodressedmuons").dressedLeptons(); _neutrinos = apply(event, "neutrinos").particlesByPt(); const Jets& all_jets = apply( event, "jets").jetsByPt(Cuts::pT > 25.0*GeV && Cuts::abseta < 2.5); //get true l+jets events by removing events with more than 1 electron||muon neutrino unsigned int n_elmu_neutrinos = 0; foreach (const Particle p, _neutrinos) { if (p.abspid() == 12 || p.abspid() == 14) ++n_elmu_neutrinos; } if (n_elmu_neutrinos != 1) vetoEvent; DressedLepton *lepton; if ( _dressedelectrons.size()) lepton = &_dressedelectrons[0]; else if (_dressedmuons.size()) lepton = &_dressedmuons[0]; else vetoEvent; // Calculate the missing ET, using the prompt neutrinos only (really?) /// @todo Why not use MissingMomentum? FourMomentum met; foreach (const Particle& p, _neutrinos) met += p.momentum(); //remove jets if they are within dR < 0.2 of lepton Jets jets; foreach(const Jet& jet, all_jets) { bool keep = true; foreach (const DressedLepton& el, _vetodressedelectrons) { keep &= deltaR(jet, el) >= 0.2; } if (keep) jets += jet; } bool overlap = false; Jets bjets, lightjets; for (unsigned int i = 0; i < jets.size(); ++i) { const Jet& jet = jets[i]; foreach (const DressedLepton& el, _dressedelectrons) overlap |= deltaR(jet, el) < 0.4; foreach (const DressedLepton& mu, _dressedmuons) overlap |= deltaR(jet, mu) < 0.4; for (unsigned int j = i + 1; j < jets.size(); ++j) { overlap |= deltaR(jet, jets[j]) < 0.5; } //// Count the number of b-tags bool b_tagged = false; // This is closer to the Particles bTags = jet.bTags(); // analysis. Something foreach ( Particle b, bTags ) { // about ghost-associated b_tagged |= b.pT() > 5*GeV; // B-hadrons } // if ( b_tagged ) bjets += jet; else lightjets += jet; } // remove events with object overlap if (overlap) vetoEvent; if (bjets.size() < 2 || lightjets.size() < 2) vetoEvent; FourMomentum pbjet1; //Momentum of bjet1 FourMomentum pbjet2; //Momentum of bjet2 if ( deltaR(bjets[0], *lepton) <= deltaR(bjets[1], *lepton) ) { pbjet1 = bjets[0].momentum(); pbjet2 = bjets[1].momentum(); } else { pbjet1 = bjets[1].momentum(); pbjet2 = bjets[0].momentum(); } FourMomentum pjet1; // Momentum of jet1 if (lightjets.size()) pjet1 = lightjets[0].momentum(); FourMomentum pjet2; // Momentum of jet 2 if (lightjets.size() > 1) pjet2 = lightjets[1].momentum(); double pz = computeneutrinoz(lepton->momentum(), met); FourMomentum ppseudoneutrino( sqrt(sqr(met.px()) + sqr(met.py()) + sqr(pz)), met.px(), met.py(), pz); //compute leptonic, hadronic, combined pseudo-top FourMomentum ppseudotoplepton = lepton->momentum() + ppseudoneutrino + pbjet1; FourMomentum ppseudotophadron = pbjet2 + pjet1 + pjet2; FourMomentum pttbar = ppseudotoplepton + ppseudotophadron; // Evaluate basic event selection bool pass_eljets = (_dressedelectrons.size() == 1) && (_vetodressedelectrons.size() < 2) && (_vetodressedmuons.empty()) && (met.pT() > 30*GeV) && (_mT(_dressedelectrons[0].momentum(), met) > 35*GeV) && (jets.size() >= 4); bool pass_mujets = (_dressedmuons.size() == 1) && (_vetodressedmuons.size() < 2) && (_vetodressedelectrons.empty()) && (met.pT() > 30*GeV) && (_mT(_dressedmuons[0].momentum(), met) > 35*GeV) && (jets.size() >= 4); // basic event selection requirements if (!pass_eljets && !pass_mujets) vetoEvent; // Fill histograms - const double weight = 1.0; - //pseudotop hadrons and leptons fill histogram - _h["ptpseudotoplepton"]->fill( ppseudotoplepton.pt(), weight); //pT of pseudo top lepton - _h["absrappseudotoplepton"]->fill(ppseudotoplepton.absrap(), weight); //absolute rapidity of pseudo top lepton - _h["ptpseudotophadron"]->fill( ppseudotophadron.pt(), weight); //pT of pseudo top hadron - _h["absrappseudotophadron"]->fill(ppseudotophadron.absrap(), weight); //absolute rapidity of pseudo top hadron - _h["absrapttbar"]->fill( pttbar.absrap(), weight); //absolute rapidity of ttbar - _h["ttbarmass"]->fill( pttbar.mass(), weight); //mass of ttbar - _h["ptttbar"]->fill( pttbar.pt(), weight); //fill pT of ttbar in combined channel + _h["ptpseudotoplepton"]->fill( ppseudotoplepton.pt()); //pT of pseudo top lepton + _h["absrappseudotoplepton"]->fill(ppseudotoplepton.absrap()); //absolute rapidity of pseudo top lepton + _h["ptpseudotophadron"]->fill( ppseudotophadron.pt()); //pT of pseudo top hadron + _h["absrappseudotophadron"]->fill(ppseudotophadron.absrap()); //absolute rapidity of pseudo top hadron + _h["absrapttbar"]->fill( pttbar.absrap()); //absolute rapidity of ttbar + _h["ttbarmass"]->fill( pttbar.mass()); //mass of ttbar + _h["ptttbar"]->fill( pttbar.pt()); //fill pT of ttbar in combined channel if (pass_eljets) { // electron channel fill histogram - _h["ptpseudotoplepton_el"]->fill( ppseudotoplepton.pt(), weight); //pT of pseudo top lepton - _h["absrappseudotoplepton_el"]->fill(ppseudotoplepton.absrap(), weight); //absolute rapidity of pseudo top lepton - _h["ptpseudotophadron_el"]->fill( ppseudotophadron.pt(), weight); //pT of pseudo top hadron - _h["absrappseudotophadron_el"]->fill(ppseudotophadron.absrap(), weight); //absolute rapidity of pseudo top hadron - _h["absrapttbar_el"]->fill( pttbar.absrap(), weight); //absolute rapidity of ttbar - _h["ttbarmass_el"]->fill( pttbar.mass(), weight); // fill electron channel ttbar mass - _h["ptttbar_el"]->fill( pttbar.pt(), weight); //fill pT of ttbar in electron channel + _h["ptpseudotoplepton_el"]->fill( ppseudotoplepton.pt()); //pT of pseudo top lepton + _h["absrappseudotoplepton_el"]->fill(ppseudotoplepton.absrap()); //absolute rapidity of pseudo top lepton + _h["ptpseudotophadron_el"]->fill( ppseudotophadron.pt()); //pT of pseudo top hadron + _h["absrappseudotophadron_el"]->fill(ppseudotophadron.absrap()); //absolute rapidity of pseudo top hadron + _h["absrapttbar_el"]->fill( pttbar.absrap()); //absolute rapidity of ttbar + _h["ttbarmass_el"]->fill( pttbar.mass()); // fill electron channel ttbar mass + _h["ptttbar_el"]->fill( pttbar.pt()); //fill pT of ttbar in electron channel } else { // muon channel fill histogram - _h["ptpseudotoplepton_mu"]->fill( ppseudotoplepton.pt(), weight); //pT of pseudo top lepton - _h["absrappseudotoplepton_mu"]->fill(ppseudotoplepton.absrap(), weight); //absolute rapidity of pseudo top lepton - _h["ptpseudotophadron_mu"]->fill( ppseudotophadron.pt(), weight); //pT of pseudo top hadron - _h["absrappseudotophadron_mu"]->fill(ppseudotophadron.absrap(), weight); //absolute rapidity of pseudo top hadron - _h["absrapttbar_mu"]->fill( pttbar.absrap(), weight); //absolute rapidity of ttbar - _h["ttbarmass_mu"]->fill( pttbar.mass(), weight); //fill muon channel histograms - _h["ptttbar_mu"]->fill( pttbar.pt(), weight); //fill pT of ttbar in electron channel + _h["ptpseudotoplepton_mu"]->fill( ppseudotoplepton.pt()); //pT of pseudo top lepton + _h["absrappseudotoplepton_mu"]->fill(ppseudotoplepton.absrap()); //absolute rapidity of pseudo top lepton + _h["ptpseudotophadron_mu"]->fill( ppseudotophadron.pt()); //pT of pseudo top hadron + _h["absrappseudotophadron_mu"]->fill(ppseudotophadron.absrap()); //absolute rapidity of pseudo top hadron + _h["absrapttbar_mu"]->fill( pttbar.absrap()); //absolute rapidity of ttbar + _h["ttbarmass_mu"]->fill( pttbar.mass()); //fill muon channel histograms + _h["ptttbar_mu"]->fill( pttbar.pt()); //fill pT of ttbar in electron channel } } void finalize() { // Normalize to cross-section const double scalefactor(crossSection() / sumOfWeights()); for (map::iterator hit = _h.begin(); hit != _h.end(); ++hit) { double sf = scalefactor; if ( (hit->first).find("_") == std::string::npos ) sf *= 0.5; scale(hit->second, sf); } } private: double computeneutrinoz(const FourMomentum& lepton, FourMomentum& met) const { //computing z component of neutrino momentum given lepton and met double pzneutrino; double m_W = 80.399; // in GeV, given in the paper double k = (( sqr( m_W ) - sqr( lepton.mass() ) ) / 2 ) + (lepton.px() * met.px() + lepton.py() * met.py()); double a = sqr ( lepton.E() )- sqr ( lepton.pz() ); double b = -2*k*lepton.pz(); double c = sqr( lepton.E() ) * sqr( met.pT() ) - sqr( k ); double discriminant = sqr(b) - 4 * a * c; double quad[2] = { (- b - sqrt(discriminant)) / (2 * a), (- b + sqrt(discriminant)) / (2 * a) }; //two possible quadratic solns if (discriminant < 0) pzneutrino = - b / (2 * a); //if the discriminant is negative else { //if the discriminant is greater than or equal to zero, take the soln with smallest absolute value double absquad[2]; for (int n=0; n<2; ++n) absquad[n] = fabs(quad[n]); if (absquad[0] < absquad[1]) pzneutrino = quad[0]; else pzneutrino = quad[1]; } if ( !std::isfinite(pzneutrino) ) std::cout << "Found non-finite value" << std::endl; return pzneutrino; } double _mT(const FourMomentum &l, FourMomentum &nu) const { return sqrt( 2 * l.pT() * nu.pT() * (1 - cos(deltaPhi(l, nu))) ); } /// @name Objects that are used by the event selection decisions vector _dressedelectrons, _vetodressedelectrons, _dressedmuons, _vetodressedmuons; Particles _neutrinos; map _h; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2015_I1345452); } diff --git a/analyses/pluginATLAS/ATLAS_2015_I1364361.cc b/analyses/pluginATLAS/ATLAS_2015_I1364361.cc --- a/analyses/pluginATLAS/ATLAS_2015_I1364361.cc +++ b/analyses/pluginATLAS/ATLAS_2015_I1364361.cc @@ -1,137 +1,136 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { /// @brief Higgs differential cross-section combination between the ATLAS measurements in the yy and 4l channels /// /// Computes Higgs transverse momentum, rapidity, jet multiplicity and leading jet pT. /// /// @author Michaela Queitsch-Maitland /// @author Dag Gillberg /// @author Florian Bernlochner /// @author Sarah Heim class ATLAS_2015_I1364361 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2015_I1364361); /// Book histograms and initialise projections before the run void init() { // All final state particles const FinalState fs; declare(fs, "FS"); // Histograms with data bins book(_h_pTH_incl ,1,1,1); book(_h_yH_incl ,2,1,1); book(_h_Njets_incl ,3,1,1); book(_h_pTj1_incl ,4,1,1); } /// Perform the per-event analysis void analyze(const Event& event) { // Get the final state particles ordered by pT const Particles& fs = apply(event, "FS").particlesByPt(); // Find a stable Higgs (mandatory) const auto higgsiter = std::find_if(fs.begin(), fs.end(), [](const Particle& p){ return p.pid() == PID::HIGGSBOSON; }); if (higgsiter == fs.end()) vetoEvent; const Particle& higgs = *higgsiter; // bool stable_higgs = false; // const Particle* higgs = 0; // for (const Particle& p : fs) { // if (p.pid() == PID::HIGGSBOSON) { // stable_higgs = true; // higgs = &p; // break; // } // } // if (!stable_higgs) { // MSG_WARNING("FATAL: No stable Higgs found in event record.\n"); // vetoEvent; // } // Loop over final state particles and fill various particle vectors Particles leptons, photons, jet_ptcls; foreach ( const Particle& ptcl, fs ) { // Do not include the Higgs in jet finding! if ( ptcl.pid() == PID::HIGGSBOSON ) continue; // Neutrinos not from hadronisation if ( ptcl.isNeutrino() && !ptcl.fromHadron() ) continue; // Electrons and muons not from hadronisation if ( ( ptcl.abspid() == PID::ELECTRON || ptcl.abspid() == PID::MUON ) && !ptcl.fromHadron() ) { leptons.push_back(ptcl); continue; } // Photons not from hadronisation if ( ptcl.abspid() == PID::PHOTON && !ptcl.fromHadron() ) { photons.push_back(ptcl); continue; } // Add particle to jet inputs jet_ptcls.push_back(ptcl); } // Match FS photons to leptons within cone R=0.1 // If they are not 'dressing' photons, add to jet particle vector foreach ( const Particle& ph, photons ) { bool fsr_photon = false; foreach ( const Particle& lep, leptons ) { if ( deltaR(ph.momentum(),lep.momentum()) < 0.1 ){ fsr_photon=true; continue; } } if ( !fsr_photon ) jet_ptcls.push_back(ph); } // Let's build the jets! By hand... const PseudoJets pjs_in = mkPseudoJets(jet_ptcls); const fastjet::JetDefinition jdef(fastjet::antikt_algorithm, 0.4); const Jets alljets = mkJets(fastjet::ClusterSequence(pjs_in, jdef).inclusive_jets()); const Jets jets = sortByPt(filterBy(alljets, Cuts::pT > 30*GeV && Cuts::absrap < 4.4)); // FastJets jet_pro(FastJets::ANTIKT, 0.4); // jet_pro.calc(jet_ptcls); // Jets jets = jet_pro.jetsByPt(Cuts::pT > 30*GeV && Cuts::absrap < 4.4); - const double weight = 1.0; - _h_pTH_incl->fill(higgs.pT(), weight); - _h_yH_incl->fill(higgs.absrap(), weight); - _h_Njets_incl->fill(jets.size() > 3 ? 3 : jets.size(), weight); - _h_pTj1_incl->fill(jets.empty() ? 0 : jets[0].pT(), weight); + _h_pTH_incl->fill(higgs.pT()); + _h_yH_incl->fill(higgs.absrap()); + _h_Njets_incl->fill(jets.size() > 3 ? 3 : jets.size()); + _h_pTj1_incl->fill(jets.empty() ? 0 : jets[0].pT()); } /// Normalise histograms etc., after the run void finalize() { const double xs = crossSectionPerEvent(); scale(_h_pTH_incl, xs); scale(_h_yH_incl, xs); scale(_h_Njets_incl, xs); scale(_h_pTj1_incl, xs); } private: Histo1DPtr _h_pTH_incl, _h_yH_incl, _h_Njets_incl, _h_pTj1_incl; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2015_I1364361); } diff --git a/analyses/pluginATLAS/ATLAS_2015_I1376945.cc b/analyses/pluginATLAS/ATLAS_2015_I1376945.cc --- a/analyses/pluginATLAS/ATLAS_2015_I1376945.cc +++ b/analyses/pluginATLAS/ATLAS_2015_I1376945.cc @@ -1,223 +1,221 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { /// @brief Colour flow in hadronic top decay at 8 TeV class ATLAS_2015_I1376945 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2015_I1376945); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { const FinalState fs; PromptFinalState promptFs(fs); promptFs.acceptTauDecays(true); promptFs.acceptMuonDecays(false); IdentifiedFinalState neutrino_fs(promptFs); neutrino_fs.acceptNeutrinos(); declare(neutrino_fs, "NEUTRINO_FS"); IdentifiedFinalState Photon(fs); Photon.acceptIdPair(PID::PHOTON); IdentifiedFinalState bare_muons_fs(promptFs); bare_muons_fs.acceptIdPair(PID::MUON); IdentifiedFinalState bare_elecs_fs(promptFs); bare_elecs_fs.acceptIdPair(PID::ELECTRON); Cut lep_cuts = (Cuts::abseta < 2.5) & (Cuts::pT > 1*MeV); DressedLeptons muons(Photon, bare_muons_fs, 0.1, lep_cuts); declare(muons, "MUONS"); DressedLeptons elecs(Photon, bare_elecs_fs, 0.1, lep_cuts); declare(elecs, "ELECS"); VetoedFinalState vfs; vfs.addVetoOnThisFinalState(muons); vfs.addVetoOnThisFinalState(elecs); vfs.addVetoOnThisFinalState(neutrino_fs); FastJets fjets(vfs, FastJets::ANTIKT, 0.4); fjets.useInvisibles(); declare(fjets, "jets"); book(h_pull_all ,4,1,1); book(h_pull_charged ,5,1,1); } /// Perform the per-event analysis void analyze(const Event& event) { - const double weight = 1.0; - /************** * JETS * **************/ const Jets& allJets = apply(event, "jets").jetsByPt(Cuts::pT > 25.0*GeV && Cuts::absrap < 2.5); const vector& all_elecs = apply(event, "ELECS").dressedLeptons(); const vector& all_muons = apply(event, "MUONS").dressedLeptons(); Jets goodJets; foreach (const Jet j, allJets) { bool keep = true; foreach (const DressedLepton el, all_elecs) keep &= deltaR(j, el) >= 0.2; if (keep) goodJets += j; } if ( goodJets.size() < 4 ) vetoEvent; /**************** * LEPTONS * ****************/ vector muons, vetoMuons; foreach (const DressedLepton mu, all_muons) { bool keep = true; foreach (const Jet j, goodJets) keep &= deltaR(j, mu) >= 0.4; if (keep && mu.pt() > 15*GeV) { vetoMuons.push_back(mu); if (mu.pt() > 25*GeV) muons.push_back(mu); } } vector elecs, vetoElecs; foreach (const DressedLepton el, all_elecs) { bool keep = true; foreach (const Jet j, goodJets) keep &= deltaR(j, el) >= 0.4; if (keep && el.pt() > 15*GeV) { vetoElecs.push_back(el); if (el.pt() > 25*GeV) elecs.push_back(el); } } if (muons.empty() && elecs.empty()) vetoEvent; bool muCandidate = !( muons.size() < 1 || vetoMuons.size() > 1 || vetoElecs.size() > 0 ); bool elCandidate = !( elecs.size() < 1 || vetoElecs.size() > 1 || vetoMuons.size() > 0 ); if (!elCandidate && !muCandidate) vetoEvent; /****************************** * ELECTRON-MUON OVERLAP * ******************************/ foreach (const DressedLepton electron, elecs) { foreach (const DressedLepton muon, muons) { double d_theta = fabs(muon.theta() - electron.theta()); double d_phi = deltaPhi(muon.phi(), electron.phi()); if (d_theta < 0.005 && d_phi < 0.005) vetoEvent; } } /**************** * NEUTRINOS * ****************/ const Particles& neutrinos = apply(event, "NEUTRINO_FS").particlesByPt(); FourMomentum metVector = FourMomentum(0.,0.,0.,0.); foreach (const Particle& n, neutrinos) { metVector += n.momentum(); } double met = metVector.pt(); if (met <= 20*GeV) vetoEvent; if ( (_mT(muCandidate? muons[0] : elecs[0], metVector) + met) <= 60. ) vetoEvent; /**************** * B-JETS * ****************/ Jets bJets, wJets; foreach(Jet j, goodJets) { bool b_tagged = false; Particles bTags = j.bTags(); foreach ( Particle b, bTags ) { b_tagged |= b.pT() > 5*GeV; } if (b_tagged) bJets += j; if (!b_tagged && j.abseta() < 2.1) wJets += j; } if ( bJets.size() < 2 || wJets.size() < 2 ) vetoEvent; double pull_angle = fabs(CalculatePullAngle(wJets[0], wJets[1], 0)); - h_pull_all->fill(pull_angle / Rivet::PI, weight); + h_pull_all->fill(pull_angle / Rivet::PI); double pull_angle_charged = fabs(CalculatePullAngle(wJets[0], wJets[1], 1)); - h_pull_charged->fill(pull_angle_charged / Rivet::PI, weight); + h_pull_charged->fill(pull_angle_charged / Rivet::PI); } Vector3 CalculatePull(Jet& jet, bool &isCharged) { Vector3 pull(0.0, 0.0, 0.0); double PT = jet.pT(); Particles& constituents = jet.particles(); Particles charged_constituents; if (isCharged) { foreach (Particle p, constituents) { if (p.threeCharge() != 0) charged_constituents += p; } constituents = charged_constituents; } // calculate axis FourMomentum axis; foreach (Particle p, constituents) axis += p.momentum(); Vector3 J(axis.rap(), axis.phi(MINUSPI_PLUSPI), 0.0); // calculate pull foreach (Particle p, constituents) { Vector3 ri = Vector3(p.rap(), p.phi(MINUSPI_PLUSPI), 0.0) - J; while (ri.y() > Rivet::PI) ri.setY(ri.y() - Rivet::TWOPI); while (ri.y() < -Rivet::PI) ri.setY(ri.y() + Rivet::TWOPI); pull.setX(pull.x() + (ri.mod() * ri.x() * p.pT()) / PT); pull.setY(pull.y() + (ri.mod() * ri.y() * p.pT()) / PT); } return pull; } double CalculatePullAngle(Jet& jet1, Jet& axisjet, bool isCharged) { Vector3 pull_vector = CalculatePull(jet1, isCharged); pull_vector = Vector3(1000.*pull_vector.x(), 1000.*pull_vector.y(), 0.); double drap = axisjet.rap() - jet1.rap(); double dphi = axisjet.phi(MINUSPI_PLUSPI) - jet1.phi(MINUSPI_PLUSPI); Vector3 j2_vector(drap, dphi, 0.0); return mapAngleMPiToPi(deltaPhi(pull_vector, j2_vector)); } double _mT(const FourMomentum &l, FourMomentum &nu) const { return sqrt( 2 * l.pT() * nu.pT() * (1 - cos(deltaPhi(l, nu))) ); } /// Normalise histograms etc., after the run void finalize() { normalize(h_pull_all); normalize(h_pull_charged); } //@} private: Histo1DPtr h_pull_all; Histo1DPtr h_pull_charged; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2015_I1376945); } diff --git a/analyses/pluginATLAS/ATLAS_2015_I1387176.cc b/analyses/pluginATLAS/ATLAS_2015_I1387176.cc --- a/analyses/pluginATLAS/ATLAS_2015_I1387176.cc +++ b/analyses/pluginATLAS/ATLAS_2015_I1387176.cc @@ -1,92 +1,91 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { /// Rivet analysis class for ATLAS_2015_I1387176 dataset class ATLAS_2015_I1387176 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2015_I1387176); /// Initialization, called once before running void init() { // Projections FastJets jets(FinalState(), FastJets::ANTIKT, 0.4); jets.useInvisibles(); declare(jets, "Jets"); // Book histograms book(_hist_EEC ,1, 1, 1); book(_hist_AEEC ,2, 1, 1); // add dummy histogram for heterogenous merging string hname = "d01-x01-y01"; const Scatter2D& ref = refData(hname); hname = "d01-x01-y02"; book(_hist_dummy ,hname, ref); } void analyze(const Event& event) { - const double evtWeight = 1.0; const Jets& jets = apply(event, "Jets").jetsByPt(Cuts::pT > 50.0*GeV && Cuts::abseta < 2.5); if (jets.size() < 2) vetoEvent; if (jets[0].pT() + jets[1].pT() < 500*GeV) vetoEvent; double sumEt = 0.0; foreach (Jet j, jets) sumEt += j.E() / cosh(j.eta()); foreach (Jet j1, jets) { double et1 = j1.E() / cosh(j1.eta()); foreach (Jet j2, jets) { double et2 = j2.E() / cosh(j2.eta()); double etWeight = et1 * et2 / ( sumEt * sumEt ); double dPhi = deltaPhi(j1, j2); double cosPhi = cos(dPhi); if (cosPhi == 1.0) cosPhi = 0.9999; - _hist_EEC->fill(cosPhi, etWeight * evtWeight); - _hist_dummy->fill(cosPhi, etWeight * evtWeight); + _hist_EEC->fill(cosPhi, etWeight); + _hist_dummy->fill(cosPhi, etWeight); } } } void finalize() { scale(_hist_dummy, crossSectionPerEvent()); normalize(_hist_EEC); vector points; size_t nBins = _hist_EEC->numBins(); for (size_t k = 0; k < nBins/2; ++k) { double x = _hist_EEC->bin(k).midpoint(); double y = _hist_EEC->bin(k).height() - _hist_EEC->bin(nBins-(k+1)).height(); double ex = _hist_EEC->bin(k).xWidth()/2; double e1 = _hist_EEC->bin(k).heightErr(); double e2 = _hist_EEC->bin(nBins-(k+1)).heightErr(); double ey = sqrt( e1 * e1 + e2 * e2 ); points.push_back(Point2D(x, y, ex, ey)); } _hist_AEEC->addPoints(points); } private: Histo1DPtr _hist_EEC; Histo1DPtr _hist_dummy; Scatter2DPtr _hist_AEEC; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2015_I1387176); } diff --git a/analyses/pluginATLAS/ATLAS_2015_I1390114.cc b/analyses/pluginATLAS/ATLAS_2015_I1390114.cc --- a/analyses/pluginATLAS/ATLAS_2015_I1390114.cc +++ b/analyses/pluginATLAS/ATLAS_2015_I1390114.cc @@ -1,170 +1,168 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { // ATLAS $tt+b(b)$ at 8~TeV class ATLAS_2015_I1390114 : public Analysis { public: DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2015_I1390114); void init() { // Eta ranges Cut eta_full = Cuts::abseta < 5.0 && Cuts::pT >= 1.0*MeV; Cut eta_lep = Cuts::abseta < 2.5; // All final state particles FinalState fs(eta_full); // Get photons to dress leptons IdentifiedFinalState photons(fs); photons.acceptIdPair(PID::PHOTON); // Projection to find the electrons IdentifiedFinalState el_id(fs); el_id.acceptIdPair(PID::ELECTRON); PromptFinalState electrons(el_id); electrons.acceptTauDecays(true); declare(electrons, "electrons"); DressedLeptons dressedelectrons(photons, electrons, 0.1, eta_lep && Cuts::pT > 25*GeV, true); declare(dressedelectrons, "dressedelectrons"); DressedLeptons ewdressedelectrons(photons, electrons, 0.1, eta_full, true); // Projection to find the muons IdentifiedFinalState mu_id(fs); mu_id.acceptIdPair(PID::MUON); PromptFinalState muons(mu_id); muons.acceptTauDecays(true); declare(muons, "muons"); DressedLeptons dressedmuons(photons, muons, 0.1, eta_lep && Cuts::pT > 25*GeV, true); declare(dressedmuons, "dressedmuons"); DressedLeptons ewdressedmuons(photons, muons, 0.1, eta_full, true); // Projection to find neutrinos and produce MET IdentifiedFinalState nu_id; nu_id.acceptNeutrinos(); PromptFinalState neutrinos(nu_id); neutrinos.acceptTauDecays(true); // Jet clustering. VetoedFinalState vfs; vfs.addVetoOnThisFinalState(ewdressedelectrons); vfs.addVetoOnThisFinalState(ewdressedmuons); vfs.addVetoOnThisFinalState(neutrinos); FastJets jets(vfs,FastJets::ANTIKT, 0.4); jets.useInvisibles(); declare(jets, "jets"); book(_histo ,1,1,1); book(_ratio, 2,1,1, true); book(_aux ,"_aux", 1, 0.5, 1.5); } void analyze(const Event& event) { - const double weight = 1.0; - // Get the selected objects, using the projections. vector electrons = apply(event, "dressedelectrons").dressedLeptons(); vector muons = apply(event, "dressedmuons").dressedLeptons(); if (electrons.empty() && muons.empty()) vetoEvent; Jets jets = apply(event, "jets").jets((Cuts::pT>20*GeV) && (Cuts::abseta < 2.5), cmpMomByPt); Jets bjets; // Check overlap of jets/leptons. foreach (Jet jet, jets) { // if dR(el,jet) < 0.4 skip the event foreach (DressedLepton el, electrons) { if (deltaR(jet, el) < 0.4) vetoEvent; } // if dR(mu,jet) < 0.4 skip the event foreach (DressedLepton mu, muons) { if (deltaR(jet, mu) < 0.4) vetoEvent; } // Count the number of b-tags // We have to check that the ghost-matched B hadrons have pT > 5 GeV // By default jet.bTags() returns all B hadrons without cuts bool btagged = jet.bTags(Cuts::pT >= 5*GeV).size(); if (btagged) bjets += jet; } // Evaluate basic event selection bool pass_1lep = (electrons.size() == 1 && muons.empty()) || (muons.size() == 1 && electrons.empty()); - if (pass_1lep && bjets.size() >= 3 && jets.size() >= 5) _histo->fill(1, weight); + if (pass_1lep && bjets.size() >= 3 && jets.size() >= 5) _histo->fill(1); if (muons.size() == 1 && electrons.size() == 1 && bjets.size() >= 3) { - if (muons[0].charge() * electrons[0].charge() < 0.0) _histo->fill(2, weight); + if (muons[0].charge() * electrons[0].charge() < 0.0) _histo->fill(2); } DressedLepton *lep1 = NULL, *lep2 = NULL; bool zveto = false; if (electrons.size() == 2 && muons.empty()) { lep1 = &electrons[0]; lep2 = &electrons[1]; } else if (muons.size() == 2 && electrons.empty()) { lep1 = &muons[0]; lep2 = &muons[1]; } else if (electrons.size() == 1 && muons.size() == 1) { lep1 = &electrons[0]; lep2 = &muons[0]; zveto = true; } if (lep1 && lep2 && lep1->charge() * lep2->charge() < 0.0 && bjets.size() >= 2) { double mass = (lep1->momentum() + lep2->momentum()).mass(); bool pass_2lep = mass > 15*GeV && (zveto || !(mass > 81*GeV && mass < 101*GeV)); if ( pass_2lep && bjets.size() >= 4) { - _histo->fill(3, weight); - _histo->fill(4, weight); + _histo->fill(3); + _histo->fill(4); } if ( pass_2lep && jets.size() >= 4) { - _aux->fill(1, weight); + _aux->fill(1); } } } void finalize() { const double sf(crossSection()/sumOfWeights()/femtobarn); scale(_histo, sf); scale(_aux, sf); // construct ratio const double n = _histo->bin(3).sumW(); const double dN = _histo->bin(3).sumW2(); const double d = _aux->bin(0).sumW(); const double dD = _aux->bin(0).sumW2(); const double r = safediv(n, d); double e = sqrt( safediv(r * (1 - r), d) ); if ( _aux->effNumEntries() != _aux->numEntries() ) { // use F. James's approximation for weighted events: e = sqrt( safediv((1 - 2 * r) * dN + r * r * dD, d * d) ); } _ratio->point(0).setY(100.0 * r, 100.0 * e); // convert into percentage } private: Histo1DPtr _histo, _aux; Scatter2DPtr _ratio; }; // Hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2015_I1390114); } diff --git a/analyses/pluginATLAS/ATLAS_2015_I1393758.cc b/analyses/pluginATLAS/ATLAS_2015_I1393758.cc --- a/analyses/pluginATLAS/ATLAS_2015_I1393758.cc +++ b/analyses/pluginATLAS/ATLAS_2015_I1393758.cc @@ -1,143 +1,141 @@ #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { class ATLAS_2015_I1393758 : public Analysis { public: - + /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2015_I1393758); - + public: void init() { addProjection(FastJets(FinalState(), FastJets::ANTIKT, 0.4), "Jets"); book(forward_kappa3 ,1, 1, 1); book(forward_kappa5 ,2, 1, 1); book(forward_kappa7 ,3, 1, 1); book(central_kappa3 ,4, 1, 1); book(central_kappa5 ,5, 1, 1); book(central_kappa7 ,6, 1, 1); book(forwardRMS_kappa3, "d07-x01-y01", true); book(forwardRMS_kappa5, "d08-x01-y01", true); book(forwardRMS_kappa7, "d09-x01-y01", true); book(centralRMS_kappa3, "d10-x01-y01", true); book(centralRMS_kappa5, "d11-x01-y01", true); book(centralRMS_kappa7, "d12-x01-y01", true); } - + /// Perform the per-event analysis void analyze(const Event& event) { - const double weight = 1.0; - Jets m_goodJets = applyProjection(event, "Jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.1); if (m_goodJets.size() < 2) vetoEvent; if (m_goodJets[0].pT() < 50*GeV) vetoEvent; if (m_goodJets[1].pT() < 50*GeV) vetoEvent; if (fabs(1.0 - m_goodJets[0].pT()/m_goodJets[1].pT()) > 0.5) vetoEvent; bool check = m_goodJets[0].abseta() < m_goodJets[1].abseta(); int pos_f = int(check); int pos_c = int(!check); double kappa3_f = CalculateJetCharge(m_goodJets[pos_f], 0.3, 0.5, 1.8); double kappa5_f = CalculateJetCharge(m_goodJets[pos_f], 0.5, 0.5, 1.2); double kappa7_f = CalculateJetCharge(m_goodJets[pos_f], 0.7, 0.5, 0.9); double pT_f = m_goodJets[pos_f].pT(); - + double kappa3_c = CalculateJetCharge(m_goodJets[pos_c], 0.3, 0.5, 1.8); double kappa5_c = CalculateJetCharge(m_goodJets[pos_c], 0.5, 0.5, 1.2); double kappa7_c = CalculateJetCharge(m_goodJets[pos_c], 0.7, 0.5, 0.9); double pT_c = m_goodJets[pos_c].pT(); - forward_kappa3->fill(pT_f, kappa3_f, weight); - forward_kappa5->fill(pT_f, kappa5_f, weight); - forward_kappa7->fill(pT_f, kappa7_f, weight); + forward_kappa3->fill(pT_f, kappa3_f); + forward_kappa5->fill(pT_f, kappa5_f); + forward_kappa7->fill(pT_f, kappa7_f); - central_kappa3->fill(pT_c, kappa3_c, weight); - central_kappa5->fill(pT_c, kappa5_c, weight); - central_kappa7->fill(pT_c, kappa7_c, weight); + central_kappa3->fill(pT_c, kappa3_c); + central_kappa5->fill(pT_c, kappa5_c); + central_kappa7->fill(pT_c, kappa7_c); } double CalculateJetCharge(Jet& jet, double kappa=0.5, double pTcut=0.5, double Qmax=1.2) { double PTkap = pow(jet.momentum().pT(),kappa); double jetcharge = 0.; foreach (const Particle& p, jet.particles()) { if (p.pT() < pTcut) continue; if (p.threeCharge()) jetcharge += (p.threeCharge()/3.)*pow(p.pT(),kappa)/PTkap; } //Overflow and underflow if (jetcharge > Qmax) jetcharge = Qmax*0.9999; if (jetcharge < -Qmax) jetcharge = -Qmax*0.9999; return jetcharge; } - + /// Normalise histograms etc., after the run void finalize() { - + if (numEvents() > 2) { for (unsigned int i = 0; i < forward_kappa3->numBins(); ++i) { double stdv_fkappa3 = forward_kappa3->bin(i).numEntries() > 1? forward_kappa3->bin(i).stdDev() : 0.0; //See Eq. 3 for the factor of two: https://web.eecs.umich.edu/~fessler/papers/files/tr/stderr.pdf double yerr_fkappa3 = safediv(sqrt(forward_kappa3->bin(i).sumW2()), 2.*forward_kappa3->bin(i).sumW()); forwardRMS_kappa3->point(i).setY(stdv_fkappa3, yerr_fkappa3); double stdv_fkappa5 = forward_kappa5->bin(i).numEntries() > 1? forward_kappa5->bin(i).stdDev() : 0.0; double yerr_fkappa5 = safediv(sqrt(forward_kappa5->bin(i).sumW2()), 2.*forward_kappa5->bin(i).sumW()); forwardRMS_kappa5->point(i).setY(stdv_fkappa5, yerr_fkappa5); double stdv_fkappa7 = forward_kappa7->bin(i).numEntries() > 1? forward_kappa7->bin(i).stdDev() : 0.0; double yerr_fkappa7 = safediv(sqrt(forward_kappa7->bin(i).sumW2()), 2.*forward_kappa7->bin(i).sumW()); forwardRMS_kappa7->point(i).setY(stdv_fkappa7, yerr_fkappa7); double stdv_ckappa3 = central_kappa3->bin(i).numEntries() > 1? central_kappa3->bin(i).stdDev() : 0.0; double yerr_ckappa3 = safediv(sqrt(central_kappa3->bin(i).sumW2()), 2.*central_kappa3->bin(i).sumW()); centralRMS_kappa3->point(i).setY(stdv_ckappa3, yerr_ckappa3); double stdv_ckappa5 = central_kappa5->bin(i).numEntries() > 1? central_kappa5->bin(i).stdDev() : 0.0; double yerr_ckappa5 = safediv(sqrt(central_kappa5->bin(i).sumW2()), 2.*central_kappa5->bin(i).sumW()); centralRMS_kappa5->point(i).setY(stdv_ckappa5, yerr_ckappa5); double stdv_ckappa7 = central_kappa7->bin(i).numEntries() > 1? central_kappa7->bin(i).stdDev() : 0.0; double yerr_ckappa7 = safediv(sqrt(central_kappa7->bin(i).sumW2()), 2.*central_kappa7->bin(i).sumW()); centralRMS_kappa7->point(i).setY(stdv_ckappa7, yerr_ckappa7); } } } private: - + Profile1DPtr forward_kappa3; Profile1DPtr forward_kappa5; Profile1DPtr forward_kappa7; Profile1DPtr central_kappa3; Profile1DPtr central_kappa5; Profile1DPtr central_kappa7; Scatter2DPtr forwardRMS_kappa3; Scatter2DPtr forwardRMS_kappa5; Scatter2DPtr forwardRMS_kappa7; Scatter2DPtr centralRMS_kappa3; Scatter2DPtr centralRMS_kappa5; Scatter2DPtr centralRMS_kappa7; - + }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2015_I1393758); } diff --git a/analyses/pluginATLAS/ATLAS_2015_I1394679.cc b/analyses/pluginATLAS/ATLAS_2015_I1394679.cc --- a/analyses/pluginATLAS/ATLAS_2015_I1394679.cc +++ b/analyses/pluginATLAS/ATLAS_2015_I1394679.cc @@ -1,203 +1,202 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { /// Differential multijet cross-section measurement in different variables class ATLAS_2015_I1394679 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2015_I1394679); /// @name Analysis methods //@ /// Book histograms and initialise projections before the run void init() { // Initialise and register projections here const FinalState fs; declare(fs, "FinalState"); FastJets fj04(fs, FastJets::ANTIKT, 0.4, JetAlg::ALL_MUONS, JetAlg::DECAY_INVISIBLES); declare(fj04, "AntiKt4jets"); // Histograms book(_h["pt1"] ,1, 1, 1); book(_h["pt2"] ,2, 1, 1); book(_h["pt3"] ,3, 1, 1); book(_h["pt4"] ,4, 1, 1); book(_h["HT"] ,5, 1, 1); book(_h["M4j"] ,6, 1, 1); // Histograms with different pt/m4j cuts for (size_t i_hist = 0; i_hist < 4; ++i_hist) { book(_h["M2jratio_"+to_str(i_hist)] , 7 + i_hist, 1, 1); book(_h["dPhiMin2j_"+to_str(i_hist)] ,11 + i_hist, 1, 1); book(_h["dPhiMin3j_"+to_str(i_hist)] ,15 + i_hist, 1, 1); book(_h["dYMin2j_"+to_str(i_hist)] ,19 + i_hist, 1, 1); book(_h["dYMin3j_"+to_str(i_hist)] ,23 + i_hist, 1, 1); book(_h["dYMax2j_"+to_str(i_hist)] ,27 + i_hist, 1, 1); for (size_t ygap = 0; ygap < 4; ++ygap) { book(_h["sumPtCent_"+to_str(ygap)+to_str(i_hist)] ,31 + i_hist + ygap * 4, 1, 1); } } } /// Perform the per-event analysis void analyze(const Event& event) { const Jets& alljetskt4 = apply(event, "AntiKt4jets").jetsByPt(Cuts::pT > 60*GeV && Cuts::absrap < 2.8); // Require 4 jets with rap < 2.8 and passing pT cuts int nJets = alljetskt4.size(); if (nJets < 4) vetoEvent; Jets leadingJetskt4 = alljetskt4; leadingJetskt4.resize(4); Jet jet1(leadingJetskt4[0]), jet2(leadingJetskt4[1]), jet3(leadingJetskt4[2]), jet4(leadingJetskt4[3]); if (jet1.pT() < 100*GeV) vetoEvent; if (jet4.pT() < 64*GeV) vetoEvent; // dR cut const double dRcut = 0.65; double drmin = 9999; for (int ijet = 0; ijet < 4; ++ijet) { for (int jjet = ijet + 1; jjet < 4; ++jjet) { double myDR = deltaR(alljetskt4[ijet], alljetskt4[jjet], RAPIDITY); if (myDR < drmin) drmin = myDR; } } if (drmin < dRcut) vetoEvent; // Variables for calculation in loops over jets FourMomentum sum_alljets; double HT = 0; // scalar sum of pt of 4 leading jets double Mjj = 99999; // minimum combined mass of 2 jets double minDphi_ij = 999, minDphi_ijk = 999; // minimum azimuthal distance btw 2 & 3 jets double maxDrap_ij = -999; // maximum rapidity distance btw 2 jets double minDrap_ij = 999, minDrap_ijk = 999; // minimum rapidity distance btw 2 & 3 jets size_t maxY_i = -1, maxY_j = -1; // Loop over 4 leading jets for (size_t ij = 0; ij< 4; ++ij) { Jet& jeti = leadingJetskt4.at(ij); sum_alljets += jeti.mom(); HT += jeti.pT(); for (size_t jj = 0; jj< 4; ++jj) { if ( ij == jj ) continue; Jet& jetj = leadingJetskt4.at(jj); const double auxDphi = fabs(deltaPhi(jeti, jetj)); minDphi_ij = std::min(auxDphi, minDphi_ij); const double auxDrap = fabs(deltaRap(jeti, jetj)); minDrap_ij = std::min(auxDrap, minDrap_ij); if (auxDrap > maxDrap_ij) { maxDrap_ij = auxDrap; maxY_i = ij; maxY_j = jj; } FourMomentum sum_twojets = jeti.mom() + jetj.mom(); Mjj = std::min(Mjj, sum_twojets.mass()); for (size_t kj = 0; kj < 4; ++kj) { if (kj == ij || kj == jj) continue; Jet& jetk = leadingJetskt4.at(kj); const double auxDphi2 = auxDphi + fabs(deltaPhi(jeti, jetk)); minDphi_ijk = std::min(auxDphi2, minDphi_ijk); const double auxDrap2 = auxDrap + fabs(deltaRap(jeti, jetk)); minDrap_ijk = std::min(auxDrap2, minDrap_ijk); } } } //end loop over 4 leading jets // Combined mass of 4 leading jets const double Mjjjj = sum_alljets.mass(); // Sum of central jets pT double sumpt_twojets_cent = 0; // Scalar sum pt of central jets, with different rapidity gaps for (size_t ijet = 0; ijet < 4; ++ijet) { if (ijet == maxY_i || ijet == maxY_j) continue; // these are the forward jets sumpt_twojets_cent += leadingJetskt4.at(ijet).pT(); } // Fill histos // Mass and pT cuts in which the analysis tables are split; values are in GeV and cuts are inclusive const double m4jcuts[4] = {500, 1000, 1500, 2000}; const double pt1cutA[4] = {100, 400, 700, 1000}; const double pt1cutB[4] = {100, 250, 400, 550}; const double rapGapCut[4] = {1, 2, 3, 4}; - const double weight = 1.0; - _h["pt1"]->fill(jet1.pt(), weight); - _h["pt2"]->fill(jet2.pt(), weight); - _h["pt3"]->fill(jet3.pt(), weight); - _h["pt4"]->fill(jet4.pt(), weight); - _h["HT"] ->fill(HT, weight); - _h["M4j"]->fill(Mjjjj, weight); + _h["pt1"]->fill(jet1.pt()); + _h["pt2"]->fill(jet2.pt()); + _h["pt3"]->fill(jet3.pt()); + _h["pt4"]->fill(jet4.pt()); + _h["HT"] ->fill(HT); + _h["M4j"]->fill(Mjjjj); for (size_t i_cut = 0; i_cut < 4; ++i_cut) { const string icutstr = to_str(i_cut); if (Mjjjj > m4jcuts[i_cut]) - _h["M2jratio_"+icutstr]->fill( Mjj/Mjjjj , weight); + _h["M2jratio_"+icutstr]->fill( Mjj/Mjjjj ); if (jet1.pT() > pt1cutA[i_cut]) { - _h["dPhiMin2j_"+icutstr]->fill(minDphi_ij , weight); - _h["dPhiMin3j_"+icutstr]->fill(minDphi_ijk, weight); - _h["dYMin2j_"+icutstr]->fill(minDrap_ij , weight); - _h["dYMin3j_"+icutstr]->fill(minDrap_ijk, weight); + _h["dPhiMin2j_"+icutstr]->fill(minDphi_ij ); + _h["dPhiMin3j_"+icutstr]->fill(minDphi_ijk); + _h["dYMin2j_"+icutstr]->fill(minDrap_ij ); + _h["dYMin3j_"+icutstr]->fill(minDrap_ijk); } if (jet1.pt() > pt1cutB[i_cut]) { - _h["dYMax2j_"+icutstr]->fill( maxDrap_ij , weight); + _h["dYMax2j_"+icutstr]->fill( maxDrap_ij ); for (size_t yy = 0; yy < 4; ++yy) { if (maxDrap_ij > rapGapCut[yy]) - _h["sumPtCent_"+to_str(yy)+icutstr]->fill(sumpt_twojets_cent, weight); + _h["sumPtCent_"+to_str(yy)+icutstr]->fill(sumpt_twojets_cent); } } } //end loop over pt/m4j cuts } /// Normalise histograms etc., after the run void finalize() { const double sf = crossSection() / sumOfWeights(); /// @todo Migrate to C++11 range-for loop for (map::iterator hit = _h.begin(); hit != _h.end(); ++hit) { scale(hit->second, sf); } } //@ private: /// @name Histograms //@{ map _h; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2015_I1394679); } diff --git a/analyses/pluginATLAS/ATLAS_2015_I1394865.cc b/analyses/pluginATLAS/ATLAS_2015_I1394865.cc --- a/analyses/pluginATLAS/ATLAS_2015_I1394865.cc +++ b/analyses/pluginATLAS/ATLAS_2015_I1394865.cc @@ -1,270 +1,268 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/WFinder.hh" #include "Rivet/Projections/LeadingParticlesFinalState.hh" #include "Rivet/Projections/UnstableFinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/MergedFinalState.hh" #include "Rivet/Projections/MissingMomentum.hh" #include "Rivet/Projections/InvMassFinalState.hh" namespace Rivet { /// Inclusive 4-lepton lineshape class ATLAS_2015_I1394865 : public Analysis { public: /// Default constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2015_I1394865); void init() { FinalState fs(Cuts::abseta < 5.0); IdentifiedFinalState photon(fs, PID::PHOTON); IdentifiedFinalState bare_EL(fs, {PID::ELECTRON, -PID::ELECTRON}); IdentifiedFinalState bare_MU(fs, {PID::MUON, -PID::MUON}); // Selection 1: ZZ-> llll selection Cut etaranges_el = Cuts::abseta < 2.5 && Cuts::pT > 7*GeV; Cut etaranges_mu = Cuts::abseta < 2.7 && Cuts::pT > 6*GeV; DressedLeptons electron_sel4l(photon, bare_EL, 0.1, etaranges_el); declare(electron_sel4l, "ELECTRON_sel4l"); DressedLeptons muon_sel4l(photon, bare_MU, 0.1, etaranges_mu); declare(muon_sel4l, "MUON_sel4l"); // Both ZZ on-shell histos book(_h_ZZ_mZZ , 1, 1, 1); book(_h_ZZ_pTZZ, 2, 1, 1); } /// Do the analysis void analyze(const Event& e) { - const double weight = 1.0; - //////////////////////////////////////////////////////////////////// // Preselection of leptons for ZZ-> llll final state //////////////////////////////////////////////////////////////////// Particles leptons_sel4l; const vector& mu_sel4l = apply(e, "MUON_sel4l").dressedLeptons(); const vector& el_sel4l = apply(e, "ELECTRON_sel4l").dressedLeptons(); const vector leptonsFS_sel4l = mu_sel4l + el_sel4l; // leptonsFS_sel4l.insert( leptonsFS_sel4l.end(), mu_sel4l.begin(), mu_sel4l.end() ); // leptonsFS_sel4l.insert( leptonsFS_sel4l.end(), el_sel4l.begin(), el_sel4l.end() ); // mu: pT > 6 GeV, eta < 2.7; ele: pT > 7 GeV, eta < 2.5 for (const DressedLepton& l : leptonsFS_sel4l) { if (l.abspid() == PID::ELECTRON) leptons_sel4l.push_back(l); // REDUNDANT: if (l.pT() > 7*GeV && l.abseta() < 2.5) else if (l.abspid() == PID::MUON) leptons_sel4l.push_back(l); // REDUNDANT: if (l.pT() > 6*GeV && l.abseta() < 2.7) } ////////////////////////////////////////////////////////////////// // Exactly two opposite charged leptons ////////////////////////////////////////////////////////////////// // Calculate total 'flavour' charge double totalcharge = 0; for (const Particle& l : leptons_sel4l) totalcharge += l.pid(); // Analyze 4 lepton events if (leptons_sel4l.size() != 4 || totalcharge != 0) vetoEvent; // Identify Z states from 4 lepton pairs Zstate Z1, Z2, Z1_alt, Z2_alt; if ( !identifyZstates(Z1, Z2, Z1_alt, Z2_alt, leptons_sel4l) ) vetoEvent; const double mZ1 = Z1.mom().mass(); const double mZ2 = Z2.mom().mass(); const double mZ1_alt = Z1_alt.mom().mass(); const double mZ2_alt = Z2_alt.mom().mass(); const double pTZ1 = Z1.mom().pT(); const double pTZ2 = Z2.mom().pT(); const double mZZ = (Z1.mom() + Z2.mom()).mass(); const double pTZZ = (Z1.mom() + Z2.mom()).pT(); // Event selections // pT(Z) > 2 GeV bool pass = pTZ1 > 2*GeV && pTZ2 > 2*GeV; if (!pass) vetoEvent; // Lepton kinematics: pT > 20, 15, 10 (8 if muon) GeV int n1 = 0, n2 = 0, n3 = 0; for (Particle& l : leptons_sel4l) { if (l.pT() > 20*GeV) ++n1; if (l.pT() > 15*GeV) ++n2; if (l.pT() > 10*GeV && l.abspid() == PID::ELECTRON) ++n3; if (l.pT() > 8*GeV && l.abspid() == PID::MUON) ++n3; } pass = pass && n1>=1 && n2>=2 && n3>=3; if (!pass) vetoEvent; // Dilepton mass: 50 < mZ1 < 120 GeV, 12 < mZ2 < 120 GeV pass = pass && mZ1 > 50*GeV && mZ1 < 120*GeV; pass = pass && mZ2 > 12*GeV && mZ2 < 120*GeV; if (!pass) vetoEvent; // Lepton separation: deltaR(l, l') > 0.1 (0.2) for same- (different-) flavor leptons for (size_t i = 0; i < leptons_sel4l.size(); ++i) { for (size_t j = i + 1; j < leptons_sel4l.size(); ++j) { const Particle& l1 = leptons_sel4l[i]; const Particle& l2 = leptons_sel4l[j]; pass = pass && deltaR(l1, l2) > (l1.abspid() == l2.abspid() ? 0.1 : 0.2); if (!pass) vetoEvent; } } // J/Psi veto: m(l+l-) > 5 GeV pass = pass && mZ1 > 5*GeV && mZ2 > 5*GeV && mZ1_alt > 5*GeV && mZ2_alt > 5*GeV; if (!pass) vetoEvent; // 80 < m4l < 1000 GeV pass = pass && mZZ > 80*GeV && mZZ < 1000*GeV; if (!pass) vetoEvent; // Fill histograms - _h_ZZ_mZZ->fill(mZZ, weight); - _h_ZZ_pTZZ->fill(pTZZ, weight); + _h_ZZ_mZZ->fill(mZZ); + _h_ZZ_pTZZ->fill(pTZZ); } /// Finalize void finalize() { const double norm = crossSection()/sumOfWeights()/femtobarn; scale(_h_ZZ_mZZ, norm); scale(_h_ZZ_pTZZ, norm); } /// Generic Z candidate struct Zstate : public ParticlePair { Zstate() { } Zstate(ParticlePair _particlepair) : ParticlePair(_particlepair) { } FourMomentum mom() const { return first.momentum() + second.momentum(); } operator FourMomentum() const { return mom(); } static bool cmppT(const Zstate& lx, const Zstate& rx) { return lx.mom().pT() < rx.mom().pT(); } }; /// @brief 4l to ZZ assignment algorithm /// /// ZZ->4l pairing /// - At least two same flavour opposite sign (SFOS) lepton pairs /// - Ambiguities in pairing are resolved following the procedure /// 1. the leading Z (Z1) is choosen as the SFOS with dilepton mass closet to Z mass /// 2. the subleading Z (Z2) is choosen as the remaining SFOS dilepton pair /// /// Z1, Z2: the selected pairing /// Z1_alt, Z2_alt: the alternative pairing (the same as Z1, Z2 in 2e2m case) bool identifyZstates(Zstate& Z1, Zstate& Z2, Zstate& Z1_alt, Zstate& Z2_alt, const Particles& leptons_sel4l) { const double ZMASS = 91.1876*GeV; bool findZZ = false; Particles part_pos_el, part_neg_el, part_pos_mu, part_neg_mu; for (const Particle& l : leptons_sel4l) { if (l.abspid() == PID::ELECTRON) { if (l.pid() < 0) part_neg_el.push_back(l); if (l.pid() > 0) part_pos_el.push_back(l); } else if (l.abspid() == PID::MUON) { if (l.pid() < 0) part_neg_mu.push_back(l); if (l.pid() > 0) part_pos_mu.push_back(l); } } // eeee/mmmm channel if ((part_neg_el.size() == 2 && part_pos_el.size() == 2) || (part_neg_mu.size() == 2 && part_pos_mu.size() == 2)) { findZZ = true; Zstate Zcand_1, Zcand_2, Zcand_3, Zcand_4; Zstate Zcand_1_tmp, Zcand_2_tmp, Zcand_3_tmp, Zcand_4_tmp; if (part_neg_el.size() == 2) { // eeee Zcand_1_tmp = Zstate( ParticlePair( part_neg_el[0], part_pos_el[0] ) ); Zcand_2_tmp = Zstate( ParticlePair( part_neg_el[0], part_pos_el[1] ) ); Zcand_3_tmp = Zstate( ParticlePair( part_neg_el[1], part_pos_el[0] ) ); Zcand_4_tmp = Zstate( ParticlePair( part_neg_el[1], part_pos_el[1] ) ); } else { // mmmm Zcand_1_tmp = Zstate( ParticlePair( part_neg_mu[0], part_pos_mu[0] ) ); Zcand_2_tmp = Zstate( ParticlePair( part_neg_mu[0], part_pos_mu[1] ) ); Zcand_3_tmp = Zstate( ParticlePair( part_neg_mu[1], part_pos_mu[0] ) ); Zcand_4_tmp = Zstate( ParticlePair( part_neg_mu[1], part_pos_mu[1] ) ); } // We can have the following pairs: (Z1 + Z4) || (Z2 + Z3) // Firstly, reorder withing each quadruplet to have // - fabs(mZ1 - ZMASS) < fabs(mZ4 - ZMASS) // - fabs(mZ2 - ZMASS) < fabs(mZ3 - ZMASS) if (fabs(Zcand_1_tmp.mom().mass() - ZMASS) < fabs(Zcand_4_tmp.mom().mass() - ZMASS)) { Zcand_1 = Zcand_1_tmp; Zcand_4 = Zcand_4_tmp; } else { Zcand_1 = Zcand_4_tmp; Zcand_4 = Zcand_1_tmp; } if (fabs(Zcand_2_tmp.mom().mass() - ZMASS) < fabs(Zcand_3_tmp.mom().mass() - ZMASS)) { Zcand_2 = Zcand_2_tmp; Zcand_3 = Zcand_3_tmp; } else { Zcand_2 = Zcand_3_tmp; Zcand_3 = Zcand_2_tmp; } // We can have the following pairs: (Z1 + Z4) || (Z2 + Z3) // Secondly, select the leading and subleading Z following // 1. the leading Z (Z1) is choosen as the SFOS with dilepton mass closet to Z mass // 2. the subleading Z (Z2) is choosen as the remaining SFOS dilepton pair if (fabs(Zcand_1.mom().mass() - ZMASS) < fabs(Zcand_2.mom().mass() - ZMASS)) { Z1 = Zcand_1; Z2 = Zcand_4; Z1_alt = Zcand_2; Z2_alt = Zcand_3; } else { Z1 = Zcand_2; Z2 = Zcand_3; Z1_alt = Zcand_1; Z2_alt = Zcand_4; } } // end of eeee/mmmm channel else if (part_neg_el.size() == 1 && part_pos_el.size() == 1 && part_neg_mu.size() == 1 && part_pos_mu.size() == 1) { // 2e2m channel findZZ = true; Zstate Zcand_1, Zcand_2; Zcand_1 = Zstate( ParticlePair( part_neg_mu[0], part_pos_mu[0] ) ); Zcand_2 = Zstate( ParticlePair( part_neg_el[0], part_pos_el[0] ) ); if (fabs(Zcand_1.mom().mass() - ZMASS) < fabs(Zcand_2.mom().mass() - ZMASS)) { Z1 = Zcand_1; Z2 = Zcand_2; } else { Z1 = Zcand_2; Z2 = Zcand_1; } Z1_alt = Z1; Z2_alt = Z2; } return findZZ; } private: Histo1DPtr _h_ZZ_pTZZ, _h_ZZ_mZZ; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2015_I1394865); }