diff --git a/analyses/pluginATLAS/ATLAS_2014_I1310835.cc b/analyses/pluginATLAS/ATLAS_2014_I1310835.cc --- a/analyses/pluginATLAS/ATLAS_2014_I1310835.cc +++ b/analyses/pluginATLAS/ATLAS_2014_I1310835.cc @@ -1,267 +1,267 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" namespace Rivet { /// @brief H(125)->ZZ->4l at 8 TeV class ATLAS_2014_I1310835 : public Analysis { public: /// Default constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2014_I1310835); void init() { const FinalState fs(Cuts::abseta < 5.0); PromptFinalState photons(Cuts::abspid == PID::PHOTON); PromptFinalState bare_el(Cuts::abspid == PID::ELECTRON); PromptFinalState bare_mu(Cuts::abspid == PID::MUON); // Selection: lepton selection Cut etaranges_el = Cuts::abseta < 2.47 && Cuts::pT > 7*GeV; DressedLeptons electron_sel4l(photons, bare_el, 0.1, etaranges_el, false); declare(electron_sel4l, "electrons"); Cut etaranges_mu = Cuts::abseta < 2.7 && Cuts::pT > 6*GeV; DressedLeptons muon_sel4l(photons, bare_mu, 0.1, etaranges_mu, false); declare(muon_sel4l, "muons"); FastJets jetpro(fs, FastJets::ANTIKT, 0.4, JetAlg::Muons::NONE, JetAlg::Invisibles::NONE); declare(jetpro, "jet"); // Book histos book(_h_pt , 1, 1, 1); book(_h_rapidity , 2, 1, 1); book(_h_m34 , 3, 1, 1); book(_h_costheta , 4, 1, 1); book(_h_njets , 5, 1, 1); book(_h_leadingjetpt, 6, 1, 1); } /// Do the analysis void analyze(const Event& e) { //////////////////////////////////////////////////////////////////// // preselection of leptons for ZZ-> llll final state //////////////////////////////////////////////////////////////////// const vector& mu_sel4l = applyProjection(e, "muons").dressedLeptons(); const vector& el_sel4l = applyProjection(e, "electrons").dressedLeptons(); vector leptonsFS_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() ); ///////////////////////////////////////////////////////////////////////////// /// H->ZZ->4l pairing ///////////////////////////////////////////////////////////////////////////// size_t el_p = 0; size_t el_n = 0; size_t mu_p = 0; size_t mu_n = 0; for (const Particle& l : leptonsFS_sel4l) { if (l.abspid() == PID::ELECTRON) { if (l.pid() < 0) ++el_n; if (l.pid() > 0) ++el_p; } else if (l.abspid() == PID::MUON) { if (l.pid() < 0) ++mu_n; if (l.pid() > 0) ++mu_p; } } bool pass_sfos = ( (el_p >=2 && el_n >=2) || (mu_p >=2 && mu_n >=2) || (el_p >=1 && el_n >=1 && mu_p >=1 && mu_n >=1) ); if (!pass_sfos) vetoEvent; Zstate Z1, Z2, Zcand; size_t n_parts = leptonsFS_sel4l.size(); size_t l1_index = 0; size_t l2_index = 0; // determine Z1 first double min_mass_diff = -1; for (size_t i = 0; i < n_parts; ++i) { for (size_t j = 0; j < n_parts; ++j) { if (i >= j) continue; if (leptonsFS_sel4l[i].pid() != -1*leptonsFS_sel4l[j].pid()) continue; //only pair SFOS leptons Zcand = Zstate( ParticlePair(leptonsFS_sel4l[i], leptonsFS_sel4l[j]) ); double mass_diff = fabs( Zcand.mom().mass() - 91.1876 ); if (min_mass_diff == -1 || mass_diff < min_mass_diff) { min_mass_diff = mass_diff; Z1 = Zcand; l1_index = i; l2_index = j; } } } //determine Z2 second min_mass_diff = -1; for (size_t i = 0; i < n_parts; ++i) { if (i == l1_index || i == l2_index) continue; for (size_t j = 0; j < n_parts; ++j) { if (j == l1_index || j == l2_index || i >= j) continue; if (leptonsFS_sel4l[i].pid() != -1*leptonsFS_sel4l[j].pid()) continue; // only pair SFOS leptons Zcand = Zstate( ParticlePair(leptonsFS_sel4l[i], leptonsFS_sel4l[j]) ); double mass_diff = fabs( Zcand.mom().mass() - 91.1876 ); if (min_mass_diff == -1 || mass_diff < min_mass_diff) { min_mass_diff = mass_diff; Z2 = Zcand; } } } Particles leptons_sel4l; leptons_sel4l.push_back(Z1.first); leptons_sel4l.push_back(Z1.second); leptons_sel4l.push_back(Z2.first); leptons_sel4l.push_back(Z2.second); //////////////////////////////////////////////////////////////////////////// // Kinematic Requirements /////////////////////////////////////////////////////////////////////////// //leading lepton pT requirement std::vector lepton_pt; for (const Particle& i : leptons_sel4l) lepton_pt.push_back(i.pT() / GeV); std::sort(lepton_pt.begin(), lepton_pt.end(), [](const double pT1, const double pT2) { return pT1 > pT2; }); if (!(lepton_pt[0] > 20*GeV && lepton_pt[1] > 15*GeV && lepton_pt[2] > 10*GeV)) vetoEvent; //invariant mass requirements if (!(inRange(Z1.mom().mass(), 50*GeV, 106*GeV) && inRange(Z2.mom().mass(), 12*GeV, 115*GeV))) vetoEvent; //lepton separation requirements for (unsigned int i = 0; i < 4; ++i) { for (unsigned int j = 0; j < 4; ++j) { if (i >= j) continue; double dR = deltaR(leptons_sel4l[i], leptons_sel4l[j]); bool sameflavor = leptons_sel4l[i].abspid() == leptons_sel4l[j].abspid(); if ( sameflavor && dR < 0.1) vetoEvent; if (!sameflavor && dR < 0.2) vetoEvent; } } // J/Psi veto requirement for (unsigned int i = 0; i < 4; ++i) { for (unsigned int j = 0; j < 4; ++j) { if (i >= j) continue; if ( leptons_sel4l[i].pid() != -1*leptons_sel4l[j].pid() ) continue; if ((leptons_sel4l[i].momentum() + leptons_sel4l[j].momentum()).mass() <= 5*GeV) vetoEvent; } } // 4-lepton invariant mass requirement double m4l = (Z1.mom() + Z2.mom()).mass(); if (!(inRange(m4l, 118*GeV, 129*GeV))) vetoEvent; //////////////////////////////////////////////////////////////////////////// // Higgs observables /////////////////////////////////////////////////////////////////////////// FourMomentum Higgs = Z1.mom() + Z2.mom(); double H4l_pt = Higgs.pt()/GeV; double H4l_rapidity = Higgs.absrap(); LorentzTransform HRF_boost; - //HRF_boost.mkFrameTransformFromBeta(Higgs.boostVector()); - HRF_boost.setBetaVec(- Higgs.boostVector()); + //HRF_boost.mkFrameTransformFromBeta(Higgs.betaVec()); + HRF_boost.setBetaVec(- Higgs.betaVec()); FourMomentum Z1_in_HRF = HRF_boost.transform( Z1.mom() ); double H4l_costheta = fabs(cos( Z1_in_HRF.theta())); double H4l_m34 = Z2.mom().mass()/GeV; //////////////////////////////////////////////////////////////////////////// // Jet observables /////////////////////////////////////////////////////////////////////////// Jets jets; for (const Jet& jet : applyProjection(e, "jet").jetsByPt(Cuts::pT > 30*GeV && Cuts::absrap < 4.4)) { bool overlaps = false; for (const Particle& lep : leptonsFS_sel4l) { if (lep.abspid() != PID::ELECTRON) continue; const double dR = deltaR(lep, jet); if (dR < 0.2) { overlaps = true; break; } } if (!overlaps) jets += jet; } size_t n_jets = jets.size(); if (n_jets > 3) n_jets = 3; std::vector jet_pt; for (const Jet& i : jets) jet_pt.push_back(i.pT()/GeV); double leading_jet_pt = n_jets? jet_pt[0] : 0.; //////////////////////////////////////////////////////////////////////////// // End of H->ZZ->llll selection: now fill histograms //////////////////////////////////////////////////////////////////////////// _h_pt->fill(H4l_pt); _h_rapidity->fill(H4l_rapidity); _h_costheta->fill(H4l_costheta); _h_m34->fill(H4l_m34); _h_njets->fill(n_jets + 1); _h_leadingjetpt->fill(leading_jet_pt); } /// 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(); } }; /// Finalize void finalize() { const double norm = crossSection()/sumOfWeights()/femtobarn; std::cout << "xsec: " << crossSection() << '\n'; std::cout << "sumw: " << sumOfWeights() << '\n'; std::cout << "femb: " << femtobarn << '\n'; std::cout << "norm: " << norm << '\n'; scale(_h_pt, norm); scale(_h_rapidity, norm); scale(_h_costheta, norm); scale(_h_m34, norm); scale(_h_njets, norm); scale(_h_leadingjetpt, norm); } private: Histo1DPtr _h_pt, _h_rapidity, _h_costheta; Histo1DPtr _h_m34, _h_njets, _h_leadingjetpt; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2014_I1310835); } diff --git a/analyses/pluginATLAS/ATLAS_2018_I1646686.cc b/analyses/pluginATLAS/ATLAS_2018_I1646686.cc --- a/analyses/pluginATLAS/ATLAS_2018_I1646686.cc +++ b/analyses/pluginATLAS/ATLAS_2018_I1646686.cc @@ -1,299 +1,299 @@ // -*- 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" #include "Rivet/Projections/PartonicTops.hh" #include "Rivet/Math/LorentzTrans.hh" namespace Rivet { class ATLAS_2018_I1646686 : public Analysis { public: /// Constructor /// @brief all-hadronic ttbar at 13 TeV DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2018_I1646686); /// Book histograms and initialise projections before the run void init() { //histogram booking book(_h["inclusive"],2,1,1); bookHistograms("t_y", 1, true); bookHistograms("t_pt", 2, true); bookHistograms("t1_pt", 3); bookHistograms("t1_y", 4); bookHistograms("t2_pt", 5); bookHistograms("t2_y", 6); bookHistograms("tt_m", 7); bookHistograms("tt_pt", 8); bookHistograms("tt_y", 9); bookHistograms("tt_chi", 10); bookHistograms("tt_yboost", 11); bookHistograms("tt_pout", 12); bookHistograms("tt_dPhi", 13); bookHistograms("tt_Ht", 14); bookHistograms("tt_cosThStar", 15); // Projections Cut dressed_lep = (Cuts::abseta < 2.5) && (Cuts::pT >= 25*GeV); Cut eta_full = (Cuts::abseta < 5.0); // 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); DressedLeptons dressedelectrons(photons, electrons, 0.1, dressed_lep); declare(dressedelectrons, "elecs"); DressedLeptons ewdressedelectrons(photons, electrons, 0.1, eta_full); // Projection to find the muons IdentifiedFinalState mu_id(fs); mu_id.acceptIdPair(PID::MUON); PromptFinalState muons(mu_id); muons.acceptTauDecays(true); DressedLeptons dressedmuons(photons, muons, 0.1, dressed_lep); declare(dressedmuons, "muons"); DressedLeptons ewdressedmuons(photons, muons, 0.1, eta_full); // Projection to find neutrinos for vetoing in jets IdentifiedFinalState nu_id(fs); nu_id.acceptNeutrinos(); PromptFinalState neutrinos(nu_id); neutrinos.acceptTauDecays(true); // Jet clustering. VetoedFinalState lvfs(fs); lvfs.addVetoOnThisFinalState(mu_id); lvfs.addVetoOnThisFinalState(nu_id); VetoedFinalState vfs; vfs.addVetoOnThisFinalState(ewdressedelectrons); vfs.addVetoOnThisFinalState(ewdressedmuons); vfs.addVetoOnThisFinalState(neutrinos); FastJets jets(vfs, FastJets::ANTIKT, 0.4); jets.useInvisibles(); declare(jets, "jets"); FastJets ljets(lvfs, FastJets::ANTIKT, 1.0); ljets.useInvisibles(); declare(ljets, "ljets" ); PartonicTops partonTops; declare(partonTops, "partonicTops"); } void analyze(const Event& event) { // Parton-level top quarks const Particles partonicTops = apply( event, "partonicTops").particlesByPt(); FourMomentum top, tbar; bool foundT = false, foundTBar = false; for (const Particle& ptop : partonicTops) { const int pid = ptop.pid(); if (pid == PID::TQUARK) { top = ptop.momentum(); foundT = true; } else if (pid == -PID::TQUARK) { tbar = ptop.momentum(); foundTBar = true; } } FourMomentum t1_parton, t2_parton, ttbar_parton; if ( foundT && foundTBar ) { t1_parton = top.pT2() > tbar.pT2() ? top : tbar; t2_parton = top.pT2() > tbar.pT2() ? tbar : top; ttbar_parton = t1_parton + t2_parton; if ( t1_parton.pT() > 500*GeV && t2_parton.pT() > 350*GeV) { const double chi_parton = calcChi(t1_parton, t2_parton); const double cosThetaStar_parton = abs(calcCosThetaStar(t1_parton, t2_parton)); const double pout_parton = abs(calcPout(t1_parton, t2_parton)); const double dPhi_parton = deltaPhi(t1_parton, t2_parton); const int randomChoice = rand() % 2; const FourMomentum& randomTopParton = (randomChoice == 0) ? t1_parton : t2_parton; fillParton("t_pt", randomTopParton.pT()/GeV); fillParton("t_y", randomTopParton.absrap()); fillParton("t1_pt", t1_parton.pT()/GeV); fillParton("t1_y", t1_parton.absrap()); fillParton("t2_pt", t2_parton.pT()/GeV); fillParton("t2_y", t2_parton.absrap()); fillParton("tt_m", ttbar_parton.mass()/GeV); fillParton("tt_pt", ttbar_parton.pT()/GeV); fillParton("tt_Ht", (t1_parton.pT() + t2_parton.pT())/GeV); fillParton("tt_y", ttbar_parton.absrap()); fillParton("tt_yboost", 0.5 * abs(t1_parton.rapidity() + t2_parton.rapidity())); fillParton("tt_chi", chi_parton); fillParton("tt_cosThStar", cosThetaStar_parton); fillParton("tt_pout", pout_parton/GeV); fillParton("tt_dPhi", dPhi_parton); } } // Get and veto on dressed leptons const vector dressedElectrons = apply(event, "elecs").dressedLeptons(); const vector dressedMuons = apply(event, "muons").dressedLeptons(); if (!dressedElectrons.empty()) vetoEvent; if (!dressedMuons.empty()) vetoEvent; // Get jets const Jets& all_jets = apply( event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5); const FastJets& ljets_fj = apply( event, "ljets"); const Jets all_ljets = ljets_fj.jetsByPt(); // Trim the large-R jets Jets trimmedJets; fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2), fastjet::SelectorPtFractionMin(0.05)); for (const Jet& jet : all_ljets) trimmedJets += ljets_fj.trimJet(jet, trimmer); trimmedJets = sortByPt(trimmedJets); // Check large-R jets Jets ljets; vector b_tagged; for (const Jet& jet : trimmedJets) { if (jet.pT() < 250 * GeV) continue; if (jet.pT() > 3000 * GeV) continue; if (jet.mass() > jet.pT()) continue; if (jet.abseta() > 2.0 ) continue; ljets += jet; b_tagged += jet.bTagged(); } if (all_jets.size() < 2) vetoEvent; if (ljets.size() < 2) vetoEvent; // Identify top and anti top, compute some event variables const FourMomentum ttbar = ljets[0].momentum() + ljets[1].momentum(); const FourMomentum t1 = ljets[0].momentum(); const FourMomentum t2 = ljets[1].momentum(); const double chi = calcChi(t1, t2); const double cosThetaStar = abs(calcCosThetaStar(t1, t2)); const double pout = abs(calcPout(t1, t2)); const double dPhi = deltaPhi(t1, t2); if ( t2.pT() < 350*GeV) vetoEvent; if ( t1.pT() < 500*GeV) vetoEvent; // b-tagging for particle done on large-R jets if (!(b_tagged[0] && b_tagged[1])) vetoEvent; // Continues with signal region cuts if ( abs(t1.mass() - 172.5 * GeV) > 50*GeV ) vetoEvent; if ( abs(t2.mass() - 172.5 * GeV) > 50*GeV ) vetoEvent; _h["inclusive"]->fill(0); fillHistograms("t1_pt", t1.pT()/GeV); fillHistograms("t1_y", t1.absrap()); fillHistograms("t2_pt", t2.pT()/GeV); fillHistograms("t2_y", t2.absrap()); fillHistograms("tt_m", ttbar.mass()/TeV); fillHistograms("tt_pt", ttbar.pT()/GeV); fillHistograms("tt_Ht", (t1.pT() + t2.pT())/GeV); fillHistograms("tt_y", ttbar.absrap()); fillHistograms("tt_yboost", 0.5 * abs(t1.rapidity() + t2.rapidity())); fillHistograms("tt_chi", chi); fillHistograms("tt_cosThStar", cosThetaStar); fillHistograms("tt_pout", pout/GeV); fillHistograms("tt_dPhi", dPhi); } void finalize() { // Normalize histograms const double sf = crossSection() / sumOfWeights(); for (auto &hist : _h) { scale(hist.second, sf); if ((hist.first.find("_norm") != string::npos) && hist.second->integral(false)>0) hist.second->normalize(1.0, false); } } double calcChi(const FourMomentum& t1, const FourMomentum& t2) { double ystar = 0.5 * (t1.rapidity()-t2.rapidity()); double chi = exp( 2 * abs(ystar)); return chi; } double calcCosThetaStar(const FourMomentum& t1, const FourMomentum& t2) { FourMomentum ttbar = t1 + t2; LorentzTransform centreOfMassTrans; ttbar.setX(0); ttbar.setY(0); - centreOfMassTrans.setBetaVec( -ttbar.boostVector() ); + centreOfMassTrans.setBetaVec( -ttbar.betaVec() ); FourMomentum t1_star = centreOfMassTrans.transform(t1); double cosThetaStar = t1_star.pz()/t1_star.p3().mod(); return cosThetaStar; } double calcPout(const FourMomentum& t1, const FourMomentum& t2) { Vector3 t1V = t1.p3(); Vector3 t2V = t2.p3(); Vector3 zUnit = Vector3(0., 0., 1.); Vector3 vPerp = zUnit.cross(t1V); double pout = vPerp.dot(t2V)/vPerp.mod(); return pout; } private: map _h; //some functions for booking, filling and scaling the histograms void fillHistograms(std::string name, double value) { _h[name]->fill(value); _h[name + "_norm"]->fill(value); } void fillParton(std::string name, double value) { _h[name + "_parton"]->fill(value); _h[name + "_parton_norm"]->fill(value); } void bookHistograms(const std::string name, unsigned int index, bool onlyParton = false) { if (!onlyParton) { book(_h[name], index, 1, 1 ); book(_h[name + "_norm"], index + 13, 1, 1 ); } book(_h[name + "_parton"], index + 82, 1, 1 ); book(_h[name + "_parton_norm"], index + 97, 1, 1 ); } }; DECLARE_RIVET_PLUGIN(ATLAS_2018_I1646686); } diff --git a/analyses/pluginCMS/CMS_2016_I1413748.cc b/analyses/pluginCMS/CMS_2016_I1413748.cc --- a/analyses/pluginCMS/CMS_2016_I1413748.cc +++ b/analyses/pluginCMS/CMS_2016_I1413748.cc @@ -1,328 +1,328 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/PartonicTops.hh" namespace Rivet { /// CMS 8 TeV dilepton channel ttbar spin correlations and polarisation analysis class CMS_2016_I1413748 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(CMS_2016_I1413748); /// Book histograms and initialise projections void init() { // Complete final state FinalState fs(-MAXDOUBLE, MAXDOUBLE, 0*GeV); // Projection for dressed electrons and muons IdentifiedFinalState photons(fs); photons.acceptIdPair(PID::PHOTON); IdentifiedFinalState el_id(fs); el_id.acceptIdPair(PID::ELECTRON); PromptFinalState electrons(el_id); declare(electrons, "Electrons"); DressedLeptons dressed_electrons(photons, electrons, 0.1); declare(dressed_electrons, "DressedElectrons"); IdentifiedFinalState mu_id(fs); mu_id.acceptIdPair(PID::MUON); PromptFinalState muons(mu_id); declare(muons, "Muons"); DressedLeptons dressed_muons(photons, muons, 0.1); declare(dressed_muons, "DressedMuons"); // Parton-level top quarks declare(PartonicTops(PartonicTops::DecayMode::E_MU, false), "LeptonicPartonTops"); // Booking of histograms // This histogram is independent of the parton-level information, and is an addition to the original analysis. // It is compared to the same data as the parton-level delta_phi histogram d02-x01-y01. book(_h_dphidressedleptons, "d00-x01-y01", _bins_dphi); // The remaining histos use parton-level information book(_h_dphi, "d02-x01-y01", _bins_dphi); book(_h_cos_opening_angle, "d05-x01-y01", _bins_cos_opening_angle); book(_h_c1c2, "d08-x01-y01", _bins_c1c2); book(_h_lep_costheta, "d11-x01-y01", _bins_lep_costheta); book(_h_lep_costheta_CPV, "d14-x01-y01", _bins_lep_costheta_CPV); // 2D histos book(_h_dphi_var[0], "d20-x01-y01", _bins_dphi, _bins_tt_mass); book(_h_cos_opening_angle_var[0], "d26-x01-y01", _bins_cos_opening_angle, _bins_tt_mass); book(_h_c1c2_var[0], "d32-x01-y01", _bins_c1c2, _bins_tt_mass); book(_h_lep_costheta_var[0], "d38-x01-y01", _bins_lep_costheta, _bins_tt_mass); book(_h_lep_costheta_CPV_var[0], "d44-x01-y01", _bins_lep_costheta_CPV, _bins_tt_mass); book(_h_dphi_var[1], "d50-x01-y01", _bins_dphi, _bins_tt_pT); book(_h_cos_opening_angle_var[1], "d56-x01-y01", _bins_cos_opening_angle, _bins_tt_pT); book(_h_c1c2_var[1], "d62-x01-y01", _bins_c1c2, _bins_tt_pT); book(_h_lep_costheta_var[1], "d68-x01-y01", _bins_lep_costheta, _bins_tt_pT); book(_h_lep_costheta_CPV_var[1], "d74-x01-y01", _bins_lep_costheta_CPV, _bins_tt_pT); book(_h_dphi_var[2], "d80-x01-y01", _bins_dphi, _bins_tt_absrapidity); book(_h_cos_opening_angle_var[2], "d86-x01-y01", _bins_cos_opening_angle, _bins_tt_absrapidity); book(_h_c1c2_var[2], "d92-x01-y01", _bins_c1c2, _bins_tt_absrapidity); book(_h_lep_costheta_var[2], "d98-x01-y01", _bins_lep_costheta, _bins_tt_absrapidity); book(_h_lep_costheta_CPV_var[2], "d104-x01-y01", _bins_lep_costheta_CPV, _bins_tt_absrapidity); // Profile histos for asymmetries book(_h_dphi_profile[0], "d17-x01-y01", _bins_tt_mass); book(_h_cos_opening_angle_profile[0], "d23-x01-y01", _bins_tt_mass); book(_h_c1c2_profile[0], "d29-x01-y01", _bins_tt_mass); book(_h_lep_costheta_profile[0], "d35-x01-y01", _bins_tt_mass); book(_h_lep_costheta_CPV_profile[0], "d41-x01-y01", _bins_tt_mass); book(_h_dphi_profile[1], "d47-x01-y01", _bins_tt_pT); book(_h_cos_opening_angle_profile[1], "d53-x01-y01", _bins_tt_pT); book(_h_c1c2_profile[1], "d59-x01-y01", _bins_tt_pT); book(_h_lep_costheta_profile[1], "d65-x01-y01", _bins_tt_pT); book(_h_lep_costheta_CPV_profile[1], "d71-x01-y01", _bins_tt_pT); book(_h_dphi_profile[2], "d77-x01-y01", _bins_tt_absrapidity); book(_h_cos_opening_angle_profile[2], "d83-x01-y01", _bins_tt_absrapidity); book(_h_c1c2_profile[2], "d89-x01-y01", _bins_tt_absrapidity); book(_h_lep_costheta_profile[2], "d95-x01-y01", _bins_tt_absrapidity); book(_h_lep_costheta_CPV_profile[2], "d101-x01-y01", _bins_tt_absrapidity); } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = 1.0; // Use particle-level leptons for the first histogram const DressedLeptons& dressed_electrons = applyProjection(event, "DressedElectrons"); const DressedLeptons& dressed_muons = applyProjection(event, "DressedMuons"); const vector dressedels = dressed_electrons.dressedLeptons(); const vector dressedmus = dressed_muons.dressedLeptons(); const size_t ndressedel = dressedels.size(); const size_t ndressedmu = dressedmus.size(); // For the particle-level histogram, require exactly one electron and exactly one muon, to select // the ttbar->emu channel. Note this means ttbar->emu events with additional PromptFinalState // dilepton pairs from the shower are vetoed - for PYTHIA8, this affects ~0.5% of events, so the // effect is well below the level of sensitivity of the measured distribution. if ( ndressedel == 1 && ndressedmu == 1 ) { const int electrontouse = 0, muontouse = 0; // Opposite-charge leptons only if ( sameSign(dressedels[electrontouse],dressedmus[muontouse]) ) { MSG_INFO("Error, e and mu have same charge, skipping event"); } else { //Get the four-momenta of the positively- and negatively-charged leptons FourMomentum lepPlus = dressedels[electrontouse].charge() > 0 ? dressedels[electrontouse] : dressedmus[muontouse]; FourMomentum lepMinus = dressedels[electrontouse].charge() > 0 ? dressedmus[muontouse] : dressedels[electrontouse]; // Now calculate the variable double dphi_temp = deltaPhi(lepPlus,lepMinus); fillWithUFOF( _h_dphidressedleptons, dphi_temp, weight ); } } // The remaining variables use parton-level information. // Get the leptonically decaying tops const Particles& leptonicpartontops = apply(event, "LeptonicPartonTops").particlesByPt(); Particles chargedleptons; unsigned int ntrueleptonictops = 0; bool oppositesign = false; if ( leptonicpartontops.size() == 2 ) { for (size_t k = 0; k < leptonicpartontops.size(); ++k) { // Get the lepton const Particle lepTop = leptonicpartontops[k]; const auto isPromptChargedLepton = [](const Particle& p){return (isChargedLepton(p) && isPrompt(p, false, false));}; Particles lepton_candidates = lepTop.allDescendants(firstParticleWith(isPromptChargedLepton), false); if ( lepton_candidates.size() < 1 ) MSG_WARNING("error, PartonicTops::DecayMode::E_MU top quark had no daughter lepton candidate, skipping event."); // In some cases there is no lepton from the W decay but only leptons from the decay of a radiated gamma. // These hadronic PartonicTops are currently being mistakenly selected by PartonicTops::DecayMode::E_MU (as of April 2017), and need to be rejected. // PartonicTops::DecayMode::E_MU is being fixed in Rivet, and when it is the veto below should do nothing. /// @todo Should no longer be necessary -- remove bool istrueleptonictop = false; for (size_t i = 0; i < lepton_candidates.size(); ++i) { const Particle& lepton_candidate = lepton_candidates[i]; if ( lepton_candidate.hasParent(PID::PHOTON) ) { MSG_DEBUG("Found gamma parent, top: " << k+1 << " of " << leptonicpartontops.size() << " , lepton: " << i+1 << " of " << lepton_candidates.size()); continue; } if ( !istrueleptonictop && sameSign(lepTop,lepton_candidate) ) { chargedleptons.push_back(lepton_candidate); istrueleptonictop = true; } else MSG_WARNING("Found extra prompt charged lepton from top decay (and without gamma parent), ignoring it."); } if ( istrueleptonictop ) ++ntrueleptonictops; } } if ( ntrueleptonictops == 2 ) { oppositesign = !( sameSign(chargedleptons[0],chargedleptons[1]) ); if ( !oppositesign ) MSG_WARNING("error, same charge tops, skipping event."); } if ( ntrueleptonictops == 2 && oppositesign ) { // Get the four-momenta of the positively- and negatively-charged leptons FourMomentum lepPlus = chargedleptons[0].charge() > 0 ? chargedleptons[0] : chargedleptons[1]; FourMomentum lepMinus = chargedleptons[0].charge() > 0 ? chargedleptons[1] : chargedleptons[0]; const double dphi_temp = deltaPhi(lepPlus,lepMinus); // Get the four-momenta of the positively- and negatively-charged tops FourMomentum topPlus_p4 = leptonicpartontops[0].pdgId() > 0 ? leptonicpartontops[0] : leptonicpartontops[1]; FourMomentum topMinus_p4 = leptonicpartontops[0].pdgId() > 0 ? leptonicpartontops[1] : leptonicpartontops[0]; const FourMomentum ttbar_p4 = topPlus_p4 + topMinus_p4; const double tt_mass_temp = ttbar_p4.mass(); const double tt_absrapidity_temp = ttbar_p4.absrapidity(); const double tt_pT_temp = ttbar_p4.pT(); // Lorentz transformations to calculate the spin observables in the helicity basis // Transform everything to the ttbar CM frame LorentzTransform ttCM; - ttCM.setBetaVec(-ttbar_p4.boostVector()); + ttCM.setBetaVec(-ttbar_p4.betaVec()); topPlus_p4 = ttCM.transform(topPlus_p4); topMinus_p4 = ttCM.transform(topMinus_p4); lepPlus = ttCM.transform(lepPlus); lepMinus = ttCM.transform(lepMinus); // Now boost the leptons to their parent top CM frames LorentzTransform topPlus, topMinus; - topPlus.setBetaVec(-topPlus_p4.boostVector()); - topMinus.setBetaVec(-topMinus_p4.boostVector()); + topPlus.setBetaVec(-topPlus_p4.betaVec()); + topMinus.setBetaVec(-topMinus_p4.betaVec()); lepPlus = topPlus.transform(lepPlus); lepMinus = topMinus.transform(lepMinus); const double lepPlus_costheta_temp = lepPlus.vector3().dot(topPlus_p4.vector3()) / (lepPlus.vector3().mod() * topPlus_p4.vector3().mod()); const double lepMinus_costheta_temp = lepMinus.vector3().dot(topMinus_p4.vector3()) / (lepMinus.vector3().mod() * topMinus_p4.vector3().mod()); const double c1c2_temp = lepPlus_costheta_temp * lepMinus_costheta_temp; const double cos_opening_angle_temp = lepPlus.vector3().dot(lepMinus.vector3()) / (lepPlus.vector3().mod() * lepMinus.vector3().mod()); // Fill parton-level histos fillWithUFOF( _h_dphi, dphi_temp, weight ); fillWithUFOF( _h_cos_opening_angle, cos_opening_angle_temp, weight ); fillWithUFOF( _h_c1c2, c1c2_temp, weight ); fillWithUFOF( _h_lep_costheta, lepPlus_costheta_temp, weight ); fillWithUFOF( _h_lep_costheta, lepMinus_costheta_temp, weight ); fillWithUFOF( _h_lep_costheta_CPV, lepPlus_costheta_temp, weight ); fillWithUFOF( _h_lep_costheta_CPV, -lepMinus_costheta_temp, weight ); // Now fill the same variables in the 2D and profile histos vs ttbar invariant mass, pT, and absolute rapidity for (int i_var = 0; i_var < 3; ++i_var) { double var; if ( i_var == 0 ) { var = tt_mass_temp; } else if ( i_var == 1 ) { var = tt_pT_temp; } else { var = tt_absrapidity_temp; } fillWithUFOF( _h_dphi_var[i_var], dphi_temp, var, weight ); fillWithUFOF( _h_cos_opening_angle_var[i_var], cos_opening_angle_temp, var, weight ); fillWithUFOF( _h_c1c2_var[i_var], c1c2_temp, var, weight ); fillWithUFOF( _h_lep_costheta_var[i_var], lepPlus_costheta_temp, var, weight ); fillWithUFOF( _h_lep_costheta_var[i_var], lepMinus_costheta_temp, var, weight ); fillWithUFOF( _h_lep_costheta_CPV_var[i_var], lepPlus_costheta_temp, var, weight ); fillWithUFOF( _h_lep_costheta_CPV_var[i_var], -lepMinus_costheta_temp, var, weight ); fillWithUFOF( _h_dphi_profile[i_var], dphi_temp, var, weight, (_h_dphi->xMax() + _h_dphi->xMin())/2. ); fillWithUFOF( _h_cos_opening_angle_profile[i_var], cos_opening_angle_temp, var, weight, (_h_cos_opening_angle->xMax() + _h_cos_opening_angle->xMin())/2. ); fillWithUFOF( _h_c1c2_profile[i_var], c1c2_temp, var, weight, (_h_c1c2->xMax() + _h_c1c2->xMin())/2. ); fillWithUFOF( _h_lep_costheta_profile[i_var], lepPlus_costheta_temp, var, weight, (_h_lep_costheta->xMax() + _h_lep_costheta->xMin())/2. ); fillWithUFOF( _h_lep_costheta_profile[i_var], lepMinus_costheta_temp, var, weight, (_h_lep_costheta->xMax() + _h_lep_costheta->xMin())/2. ); fillWithUFOF( _h_lep_costheta_CPV_profile[i_var], lepPlus_costheta_temp, var, weight, (_h_lep_costheta_CPV->xMax() + _h_lep_costheta_CPV->xMin())/2. ); fillWithUFOF( _h_lep_costheta_CPV_profile[i_var], -lepMinus_costheta_temp, var, weight, (_h_lep_costheta_CPV->xMax() + _h_lep_costheta_CPV->xMin())/2. ); } } } /// Normalise histograms to unit area void finalize() { normalize(_h_dphidressedleptons); normalize(_h_dphi); normalize(_h_cos_opening_angle); normalize(_h_c1c2); normalize(_h_lep_costheta); normalize(_h_lep_costheta_CPV); for (int i_var = 0; i_var < 3; ++i_var) { normalize(_h_dphi_var[i_var]); normalize(_h_cos_opening_angle_var[i_var]); normalize(_h_c1c2_var[i_var]); normalize(_h_lep_costheta_var[i_var]); normalize(_h_lep_costheta_CPV_var[i_var]); } } private: Histo1DPtr _h_dphidressedleptons, _h_dphi, _h_lep_costheta, _h_lep_costheta_CPV, _h_c1c2, _h_cos_opening_angle; Histo2DPtr _h_dphi_var[3], _h_lep_costheta_var[3], _h_lep_costheta_CPV_var[3], _h_c1c2_var[3], _h_cos_opening_angle_var[3]; Profile1DPtr _h_dphi_profile[3], _h_lep_costheta_profile[3], _h_lep_costheta_CPV_profile[3], _h_c1c2_profile[3], _h_cos_opening_angle_profile[3]; const vector _bins_tt_mass = {300., 430., 530., 1200.}; const vector _bins_tt_pT = {0., 41., 92., 300.}; const vector _bins_tt_absrapidity = {0., 0.34, 0.75, 1.5}; const vector _bins_dphi = {0., 5.*M_PI/60., 10.*M_PI/60., 15.*M_PI/60., 20.*M_PI/60., 25.*M_PI/60., 30.*M_PI/60., 35.*M_PI/60., 40.*M_PI/60., 45.*M_PI/60., 50.*M_PI/60., 55.*M_PI/60., M_PI}; const vector _bins_lep_costheta = {-1., -2./3., -1./3., 0., 1./3., 2./3., 1.}; const vector _bins_lep_costheta_CPV = {-1., -2./3., -1./3., 0., 1./3., 2./3., 1.}; const vector _bins_c1c2 = {-1., -0.4, -10./60., 0., 10./60., 0.4, 1.}; const vector _bins_cos_opening_angle = {-1., -2./3., -1./3., 0., 1./3., 2./3., 1.}; void fillWithUFOF(Histo1DPtr h, double x, double w) { h->fill(std::max(std::min(x, h->xMax()-1e-9),h->xMin()+1e-9), w); } void fillWithUFOF(Histo2DPtr h, double x, double y, double w) { h->fill(std::max(std::min(x, h->xMax()-1e-9),h->xMin()+1e-9), std::max(std::min(y, h->yMax()-1e-9),h->yMin()+1e-9), w); } void fillWithUFOF(Profile1DPtr h, double x, double y, double w, double c) { h->fill(std::max(std::min(y, h->xMax()-1e-9),h->xMin()+1e-9), float(x > c) - float(x < c), w); } }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_2016_I1413748); } diff --git a/include/Rivet/Math/Vector4.hh b/include/Rivet/Math/Vector4.hh --- a/include/Rivet/Math/Vector4.hh +++ b/include/Rivet/Math/Vector4.hh @@ -1,1540 +1,1536 @@ #ifndef RIVET_MATH_VECTOR4 #define RIVET_MATH_VECTOR4 #include "Rivet/Math/MathHeader.hh" #include "Rivet/Math/MathUtils.hh" #include "Rivet/Math/VectorN.hh" #include "Rivet/Math/Vector3.hh" namespace Rivet { class FourVector; class FourMomentum; class LorentzTransform; typedef FourVector Vector4; FourVector transform(const LorentzTransform& lt, const FourVector& v4); /// @brief Specialisation of VectorN to a general (non-momentum) Lorentz 4-vector. /// /// @todo Add composite set/mk methods from different coord systems class FourVector : public Vector<4> { friend FourVector multiply(const double a, const FourVector& v); friend FourVector multiply(const FourVector& v, const double a); friend FourVector add(const FourVector& a, const FourVector& b); friend FourVector transform(const LorentzTransform& lt, const FourVector& v4); public: FourVector() : Vector<4>() { } template FourVector(const V4& other) { this->setT(other.t()); this->setX(other.x()); this->setY(other.y()); this->setZ(other.z()); } FourVector(const Vector<4>& other) : Vector<4>(other) { } FourVector(const double t, const double x, const double y, const double z) { this->setT(t); this->setX(x); this->setY(y); this->setZ(z); } virtual ~FourVector() { } public: double t() const { return get(0); } double t2() const { return sqr(t()); } FourVector& setT(const double t) { set(0, t); return *this; } double x() const { return get(1); } double x2() const { return sqr(x()); } FourVector& setX(const double x) { set(1, x); return *this; } double y() const { return get(2); } double y2() const { return sqr(y()); } FourVector& setY(const double y) { set(2, y); return *this; } double z() const { return get(3); } double z2() const { return sqr(z()); } FourVector& setZ(const double z) { set(3, z); return *this; } double invariant() const { // Done this way for numerical precision return (t() + z())*(t() - z()) - x()*x() - y()*y(); } bool isNull() const { return Rivet::isZero(invariant()); } /// Angle between this vector and another double angle(const FourVector& v) const { return vector3().angle( v.vector3() ); } /// Angle between this vector and another (3-vector) double angle(const Vector3& v3) const { return vector3().angle(v3); } /// @brief Mod-square of the projection of the 3-vector on to the \f$ x-y \f$ plane /// This is a more efficient function than @c polarRadius, as it avoids the square root. /// Use it if you only need the squared value, or e.g. an ordering by magnitude. double polarRadius2() const { return vector3().polarRadius2(); } /// Synonym for polarRadius2 double perp2() const { return vector3().perp2(); } /// Synonym for polarRadius2 double rho2() const { return vector3().rho2(); } /// Magnitude of projection of 3-vector on to the \f$ x-y \f$ plane double polarRadius() const { return vector3().polarRadius(); } /// Synonym for polarRadius double perp() const { return vector3().perp(); } /// Synonym for polarRadius double rho() const { return vector3().rho(); } /// Projection of 3-vector on to the \f$ x-y \f$ plane Vector3 polarVec() const { return vector3().polarVec(); } /// Synonym for polarVec Vector3 perpVec() const { return vector3().perpVec(); } /// Synonym for polarVec Vector3 rhoVec() const { return vector3().rhoVec(); } /// Angle subtended by the 3-vector's projection in x-y and the x-axis. double azimuthalAngle(const PhiMapping mapping=ZERO_2PI) const { return vector3().azimuthalAngle(mapping); } /// Synonym for azimuthalAngle. double phi(const PhiMapping mapping=ZERO_2PI) const { return vector3().phi(mapping); } /// Angle subtended by the 3-vector and the z-axis. double polarAngle() const { return vector3().polarAngle(); } /// Synonym for polarAngle. double theta() const { return vector3().theta(); } /// Pseudorapidity (defined purely by the 3-vector components) double pseudorapidity() const { return vector3().pseudorapidity(); } /// Synonym for pseudorapidity. double eta() const { return vector3().eta(); } /// Get the \f$ |\eta| \f$ directly. double abspseudorapidity() const { return fabs(eta()); } /// Get the \f$ |\eta| \f$ directly (alias). double abseta() const { return fabs(eta()); } /// Get the spatial part of the 4-vector as a 3-vector. Vector3 vector3() const { return Vector3(get(1), get(2), get(3)); } /// Implicit cast to a 3-vector operator Vector3 () const { return vector3(); } public: /// Contract two 4-vectors, with metric signature (+ - - -). double contract(const FourVector& v) const { const double result = t()*v.t() - x()*v.x() - y()*v.y() - z()*v.z(); return result; } /// Contract two 4-vectors, with metric signature (+ - - -). double dot(const FourVector& v) const { return contract(v); } /// Contract two 4-vectors, with metric signature (+ - - -). double operator*(const FourVector& v) const { return contract(v); } /// Multiply by a scalar. FourVector& operator*=(double a) { _vec = multiply(a, *this)._vec; return *this; } /// Divide by a scalar. FourVector& operator/=(double a) { _vec = multiply(1.0/a, *this)._vec; return *this; } /// Add to this 4-vector. FourVector& operator+=(const FourVector& v) { _vec = add(*this, v)._vec; return *this; } /// Subtract from this 4-vector. NB time as well as space components are subtracted. FourVector& operator-=(const FourVector& v) { _vec = add(*this, -v)._vec; return *this; } /// Multiply all components (space and time) by -1. FourVector operator-() const { FourVector result; result._vec = -_vec; return result; } /// Multiply space components only by -1. FourVector reverse() const { FourVector result = -*this; result.setT(-result.t()); return result; } }; /// Contract two 4-vectors, with metric signature (+ - - -). inline double contract(const FourVector& a, const FourVector& b) { return a.contract(b); } /// Contract two 4-vectors, with metric signature (+ - - -). inline double dot(const FourVector& a, const FourVector& b) { return contract(a, b); } inline FourVector multiply(const double a, const FourVector& v) { FourVector result; result._vec = a * v._vec; return result; } inline FourVector multiply(const FourVector& v, const double a) { return multiply(a, v); } inline FourVector operator*(const double a, const FourVector& v) { return multiply(a, v); } inline FourVector operator*(const FourVector& v, const double a) { return multiply(a, v); } inline FourVector operator/(const FourVector& v, const double a) { return multiply(1.0/a, v); } inline FourVector add(const FourVector& a, const FourVector& b) { FourVector result; result._vec = a._vec + b._vec; return result; } inline FourVector operator+(const FourVector& a, const FourVector& b) { return add(a, b); } inline FourVector operator-(const FourVector& a, const FourVector& b) { return add(a, -b); } /// Calculate the Lorentz self-invariant of a 4-vector. /// \f$ v_\mu v^\mu = g_{\mu\nu} x^\mu x^\nu \f$. inline double invariant(const FourVector& lv) { return lv.invariant(); } /// Angle (in radians) between spatial parts of two Lorentz vectors. inline double angle(const FourVector& a, const FourVector& b) { return a.angle(b); } /// Angle (in radians) between spatial parts of two Lorentz vectors. inline double angle(const Vector3& a, const FourVector& b) { return angle( a, b.vector3() ); } /// Angle (in radians) between spatial parts of two Lorentz vectors. inline double angle(const FourVector& a, const Vector3& b) { return a.angle(b); } //////////////////////////////////////////////// /// Specialized version of the FourVector with momentum/energy functionality. class FourMomentum : public FourVector { friend FourMomentum multiply(const double a, const FourMomentum& v); friend FourMomentum multiply(const FourMomentum& v, const double a); friend FourMomentum add(const FourMomentum& a, const FourMomentum& b); friend FourMomentum transform(const LorentzTransform& lt, const FourMomentum& v4); public: FourMomentum() { } template FourMomentum(const V4& other) { this->setE(other.t()); this->setPx(other.x()); this->setPy(other.y()); this->setPz(other.z()); } FourMomentum(const Vector<4>& other) : FourVector(other) { } FourMomentum(const double E, const double px, const double py, const double pz) { this->setE(E); this->setPx(px); this->setPy(py); this->setPz(pz); } ~FourMomentum() {} public: /// @name Coordinate setters //@{ /// Set energy \f$ E \f$ (time component of momentum). FourMomentum& setE(double E) { setT(E); return *this; } /// Set x-component of momentum \f$ p_x \f$. FourMomentum& setPx(double px) { setX(px); return *this; } /// Set y-component of momentum \f$ p_y \f$. FourMomentum& setPy(double py) { setY(py); return *this; } /// Set z-component of momentum \f$ p_z \f$. FourMomentum& setPz(double pz) { setZ(pz); return *this; } /// Set the p coordinates and energy simultaneously FourMomentum& setPE(double px, double py, double pz, double E) { if (E < 0) throw std::invalid_argument("Negative energy given as argument: " + to_str(E)); setPx(px); setPy(py); setPz(pz); setE(E); return *this; } /// Alias for setPE FourMomentum& setXYZE(double px, double py, double pz, double E) { return setPE(px, py, pz, E); } // /// Near-alias with switched arg order // FourMomentum& setEP(double E, double px, double py, double pz) { // return setPE(px, py, pz, E); // } // /// Alias for setEP // FourMomentum& setEXYZ(double E, double px, double py, double pz) { // return setEP(E, px, py, pz); // } /// Set the p coordinates and mass simultaneously FourMomentum& setPM(double px, double py, double pz, double mass) { if (mass < 0) throw std::invalid_argument("Negative mass given as argument: " + to_str(mass)); const double E = sqrt( sqr(mass) + sqr(px) + sqr(py) + sqr(pz) ); // setPx(px); setPy(py); setPz(pz); setE(E); return setPE(px, py, pz, E); } /// Alias for setPM FourMomentum& setXYZM(double px, double py, double pz, double mass) { return setPM(px, py, pz, mass); } /// Set the vector state from (eta,phi,energy) coordinates and the mass /// /// eta = -ln(tan(theta/2)) /// -> theta = 2 atan(exp(-eta)) FourMomentum& setEtaPhiME(double eta, double phi, double mass, double E) { if (mass < 0) throw std::invalid_argument("Negative mass given as argument"); if (E < 0) throw std::invalid_argument("Negative energy given as argument"); const double theta = 2 * atan(exp(-eta)); if (theta < 0 || theta > M_PI) throw std::domain_error("Polar angle outside 0..pi in calculation"); setThetaPhiME(theta, phi, mass, E); return *this; } /// Set the vector state from (eta,phi,pT) coordinates and the mass /// /// eta = -ln(tan(theta/2)) /// -> theta = 2 atan(exp(-eta)) FourMomentum& setEtaPhiMPt(double eta, double phi, double mass, double pt) { if (mass < 0) throw std::invalid_argument("Negative mass given as argument"); if (pt < 0) throw std::invalid_argument("Negative transverse momentum given as argument"); const double theta = 2 * atan(exp(-eta)); if (theta < 0 || theta > M_PI) throw std::domain_error("Polar angle outside 0..pi in calculation"); const double p = pt / sin(theta); const double E = sqrt( sqr(p) + sqr(mass) ); setThetaPhiME(theta, phi, mass, E); return *this; } /// Set the vector state from (y,phi,energy) coordinates and the mass /// /// y = 0.5 * ln((E+pz)/(E-pz)) /// -> (E^2 - pz^2) exp(2y) = (E+pz)^2 /// & (E^2 - pz^2) exp(-2y) = (E-pz)^2 /// -> E = sqrt(pt^2 + m^2) cosh(y) /// -> pz = sqrt(pt^2 + m^2) sinh(y) /// -> sqrt(pt^2 + m^2) = E / cosh(y) FourMomentum& setRapPhiME(double y, double phi, double mass, double E) { if (mass < 0) throw std::invalid_argument("Negative mass given as argument"); if (E < 0) throw std::invalid_argument("Negative energy given as argument"); const double sqrt_pt2_m2 = E / cosh(y); const double pt = sqrt( sqr(sqrt_pt2_m2) - sqr(mass) ); if (pt < 0) throw std::domain_error("Negative transverse momentum in calculation"); const double pz = sqrt_pt2_m2 * sinh(y); const double px = pt * cos(phi); const double py = pt * sin(phi); setPE(px, py, pz, E); return *this; } /// Set the vector state from (y,phi,pT) coordinates and the mass /// /// y = 0.5 * ln((E+pz)/(E-pz)) /// -> E = sqrt(pt^2 + m^2) cosh(y) [see above] FourMomentum& setRapPhiMPt(double y, double phi, double mass, double pt) { if (mass < 0) throw std::invalid_argument("Negative mass given as argument"); if (pt < 0) throw std::invalid_argument("Negative transverse mass given as argument"); const double E = sqrt( sqr(pt) + sqr(mass) ) * cosh(y); if (E < 0) throw std::domain_error("Negative energy in calculation"); setRapPhiME(y, phi, mass, E); return *this; } /// Set the vector state from (theta,phi,energy) coordinates and the mass /// /// p = sqrt(E^2 - mass^2) /// pz = p cos(theta) /// pt = p sin(theta) FourMomentum& setThetaPhiME(double theta, double phi, double mass, double E) { if (theta < 0 || theta > M_PI) throw std::invalid_argument("Polar angle outside 0..pi given as argument"); if (mass < 0) throw std::invalid_argument("Negative mass given as argument"); if (E < 0) throw std::invalid_argument("Negative energy given as argument"); const double p = sqrt( sqr(E) - sqr(mass) ); const double pz = p * cos(theta); const double pt = p * sin(theta); if (pt < 0) throw std::invalid_argument("Negative transverse momentum in calculation"); const double px = pt * cos(phi); const double py = pt * sin(phi); setPE(px, py, pz, E); return *this; } /// Set the vector state from (theta,phi,pT) coordinates and the mass /// /// p = pt / sin(theta) /// pz = p cos(theta) /// E = sqrt(p^2 + mass^2) FourMomentum& setThetaPhiMPt(double theta, double phi, double mass, double pt) { if (theta < 0 || theta > M_PI) throw std::invalid_argument("Polar angle outside 0..pi given as argument"); if (mass < 0) throw std::invalid_argument("Negative mass given as argument"); if (pt < 0) throw std::invalid_argument("Negative transverse momentum given as argument"); const double p = pt / sin(theta); const double px = pt * cos(phi); const double py = pt * sin(phi); const double pz = p * cos(theta); const double E = sqrt( sqr(p) + sqr(mass) ); setPE(px, py, pz, E); return *this; } /// Set the vector state from (pT,phi,energy) coordinates and the mass /// /// pz = sqrt(E^2 - mass^2 - pt^2) FourMomentum& setPtPhiME(double pt, double phi, double mass, double E) { if (pt < 0) throw std::invalid_argument("Negative transverse momentum given as argument"); if (mass < 0) throw std::invalid_argument("Negative mass given as argument"); if (E < 0) throw std::invalid_argument("Negative energy given as argument"); const double px = pt * cos(phi); const double py = pt * sin(phi); const double pz = sqrt(sqr(E) - sqr(mass) - sqr(pt)); setPE(px, py, pz, E); return *this; } //@} /// @name Accessors //@{ /// Get energy \f$ E \f$ (time component of momentum). double E() const { return t(); } /// Get energy-squared \f$ E^2 \f$. double E2() const { return t2(); } /// Get x-component of momentum \f$ p_x \f$. double px() const { return x(); } /// Get x-squared \f$ p_x^2 \f$. double px2() const { return x2(); } /// Get y-component of momentum \f$ p_y \f$. double py() const { return y(); } /// Get y-squared \f$ p_y^2 \f$. double py2() const { return y2(); } /// Get z-component of momentum \f$ p_z \f$. double pz() const { return z(); } /// Get z-squared \f$ p_z^2 \f$. double pz2() const { return z2(); } /// @brief Get the mass \f$ m = \sqrt{E^2 - p^2} \f$ (the Lorentz self-invariant). /// /// For spacelike momenta, the mass will be -sqrt(|mass2|). double mass() const { // assert(Rivet::isZero(mass2()) || mass2() > 0); // if (Rivet::isZero(mass2())) { // return 0.0; // } else { // return sqrt(mass2()); // } return sign(mass2()) * sqrt(fabs(mass2())); } /// Get the squared mass \f$ m^2 = E^2 - p^2 \f$ (the Lorentz self-invariant). double mass2() const { return invariant(); } /// Get 3-momentum part, \f$ p \f$. Vector3 p3() const { return vector3(); } /// Get the modulus of the 3-momentum double p() const { return p3().mod(); } /// Get the modulus-squared of the 3-momentum double p2() const { return p3().mod2(); } /// Calculate the rapidity. double rapidity() const { return 0.5 * std::log( (E() + pz()) / (E() - pz()) ); } /// Alias for rapidity. double rap() const { return rapidity(); } /// Absolute rapidity. double absrapidity() const { return fabs(rapidity()); } /// Absolute rapidity. double absrap() const { return fabs(rap()); } /// Calculate the transverse momentum vector \f$ \vec{p}_T \f$. Vector3 pTvec() const { return p3().polarVec(); } /// Synonym for pTvec Vector3 ptvec() const { return pTvec(); } /// Calculate the squared transverse momentum \f$ p_T^2 \f$. double pT2() const { return vector3().polarRadius2(); } /// Calculate the squared transverse momentum \f$ p_T^2 \f$. double pt2() const { return vector3().polarRadius2(); } /// Calculate the transverse momentum \f$ p_T \f$. double pT() const { return sqrt(pT2()); } /// Calculate the transverse momentum \f$ p_T \f$. double pt() const { return sqrt(pT2()); } /// Calculate the transverse energy \f$ E_T^2 = E^2 \sin^2{\theta} \f$. double Et2() const { return Et() * Et(); } /// Calculate the transverse energy \f$ E_T = E \sin{\theta} \f$. double Et() const { return E() * sin(polarAngle()); } //@} /// @name Lorentz boost factors and vectors //@{ /// Calculate the boost factor \f$ \gamma \f$. /// @note \f$ \gamma = E/mc^2 \f$ so we rely on the c=1 convention double gamma() const { return sqrt(E2()/mass2()); } /// Calculate the boost vector \f$ \vec{\gamma} \f$. /// @note \f$ \gamma = E/mc^2 \f$ so we rely on the c=1 convention Vector3 gammaVec() const { return gamma() * p3().unit(); } /// Calculate the boost factor \f$ \beta \f$. /// @note \f$ \beta = pc/E \f$ so we rely on the c=1 convention double beta() const { return p()/E(); } /// Calculate the boost vector \f$ \vec{\beta} \f$. /// @note \f$ \beta = pc/E \f$ so we rely on the c=1 convention Vector3 betaVec() const { // return Vector3(px()/E(), py()/E(), pz()/E()); return p3()/E(); } - /// @brief Deprecated alias for betaVec - /// @deprecated This will be removed; use betaVec() instead - Vector3 boostVector() const { return betaVec(); } - //@} //////////////////////////////////////// /// @name Sorting helpers //@{ /// Struct for sorting by increasing energy struct byEAscending { bool operator()(const FourMomentum& left, const FourMomentum& right) const{ const double pt2left = left.E(); const double pt2right = right.E(); return pt2left < pt2right; } bool operator()(const FourMomentum* left, const FourMomentum* right) const{ return (*this)(*left, *right); } }; /// Struct for sorting by decreasing energy struct byEDescending { bool operator()(const FourMomentum& left, const FourMomentum& right) const{ return byEAscending()(right, left); } bool operator()(const FourMomentum* left, const FourVector* right) const{ return (*this)(*left, *right); } }; //@} //////////////////////////////////////// /// @name Arithmetic operators //@{ /// Multiply by a scalar FourMomentum& operator*=(double a) { _vec = multiply(a, *this)._vec; return *this; } /// Divide by a scalar FourMomentum& operator/=(double a) { _vec = multiply(1.0/a, *this)._vec; return *this; } /// Add to this 4-vector. NB time as well as space components are added. FourMomentum& operator+=(const FourMomentum& v) { _vec = add(*this, v)._vec; return *this; } /// Subtract from this 4-vector. NB time as well as space components are subtracted. FourMomentum& operator-=(const FourMomentum& v) { _vec = add(*this, -v)._vec; return *this; } /// Multiply all components (time and space) by -1. FourMomentum operator-() const { FourMomentum result; result._vec = -_vec; return result; } /// Multiply space components only by -1. FourMomentum reverse() const { FourMomentum result = -*this; result.setE(-result.E()); return result; } //@} //////////////////////////////////////// /// @name Factory functions //@{ /// Make a vector from (px,py,pz,E) coordinates static FourMomentum mkXYZE(double px, double py, double pz, double E) { return FourMomentum().setPE(px, py, pz, E); } /// Make a vector from (px,py,pz) coordinates and the mass static FourMomentum mkXYZM(double px, double py, double pz, double mass) { return FourMomentum().setPM(px, py, pz, mass); } /// Make a vector from (eta,phi,energy) coordinates and the mass static FourMomentum mkEtaPhiME(double eta, double phi, double mass, double E) { return FourMomentum().setEtaPhiME(eta, phi, mass, E); } /// Make a vector from (eta,phi,pT) coordinates and the mass static FourMomentum mkEtaPhiMPt(double eta, double phi, double mass, double pt) { return FourMomentum().setEtaPhiMPt(eta, phi, mass, pt); } /// Make a vector from (y,phi,energy) coordinates and the mass static FourMomentum mkRapPhiME(double y, double phi, double mass, double E) { return FourMomentum().setRapPhiME(y, phi, mass, E); } /// Make a vector from (y,phi,pT) coordinates and the mass static FourMomentum mkRapPhiMPt(double y, double phi, double mass, double pt) { return FourMomentum().setRapPhiMPt(y, phi, mass, pt); } /// Make a vector from (theta,phi,energy) coordinates and the mass static FourMomentum mkThetaPhiME(double theta, double phi, double mass, double E) { return FourMomentum().setThetaPhiME(theta, phi, mass, E); } /// Make a vector from (theta,phi,pT) coordinates and the mass static FourMomentum mkThetaPhiMPt(double theta, double phi, double mass, double pt) { return FourMomentum().setThetaPhiMPt(theta, phi, mass, pt); } /// Make a vector from (pT,phi,energy) coordinates and the mass static FourMomentum mkPtPhiME(double pt, double phi, double mass, double E) { return FourMomentum().setPtPhiME(pt, phi, mass, E); } //@} }; inline FourMomentum multiply(const double a, const FourMomentum& v) { FourMomentum result; result._vec = a * v._vec; return result; } inline FourMomentum multiply(const FourMomentum& v, const double a) { return multiply(a, v); } inline FourMomentum operator*(const double a, const FourMomentum& v) { return multiply(a, v); } inline FourMomentum operator*(const FourMomentum& v, const double a) { return multiply(a, v); } inline FourMomentum operator/(const FourMomentum& v, const double a) { return multiply(1.0/a, v); } inline FourMomentum add(const FourMomentum& a, const FourMomentum& b) { FourMomentum result; result._vec = a._vec + b._vec; return result; } inline FourMomentum operator+(const FourMomentum& a, const FourMomentum& b) { return add(a, b); } inline FourMomentum operator-(const FourMomentum& a, const FourMomentum& b) { return add(a, -b); } ////////////////////////////////////////////////////// /// @name \f$ \Delta R \f$ calculations from 4-vectors //@{ /// @brief Calculate the squared 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors. /// /// There is a scheme ambiguity for momentum-type four vectors as to whether /// the pseudorapidity (a purely geometric concept) or the rapidity (a /// relativistic energy-momentum quantity) is to be used: this can be chosen /// via the optional scheme parameter. Use of this scheme option is /// discouraged in this case since @c RAPIDITY is only a valid option for /// vectors whose type is really the FourMomentum derived class. inline double deltaR2(const FourVector& a, const FourVector& b, RapScheme scheme=PSEUDORAPIDITY) { switch (scheme) { case PSEUDORAPIDITY : return deltaR2(a.vector3(), b.vector3()); case RAPIDITY: { const FourMomentum* ma = dynamic_cast(&a); const FourMomentum* mb = dynamic_cast(&b); if (!ma || !mb) { string err = "deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors"; throw std::runtime_error(err); } return deltaR2(*ma, *mb, scheme); } default: throw std::runtime_error("The specified deltaR scheme is not yet implemented"); } } /// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors. /// /// There is a scheme ambiguity for momentum-type four vectors as to whether /// the pseudorapidity (a purely geometric concept) or the rapidity (a /// relativistic energy-momentum quantity) is to be used: this can be chosen /// via the optional scheme parameter. Use of this scheme option is /// discouraged in this case since @c RAPIDITY is only a valid option for /// vectors whose type is really the FourMomentum derived class. inline double deltaR(const FourVector& a, const FourVector& b, RapScheme scheme=PSEUDORAPIDITY) { return sqrt(deltaR2(a, b, scheme)); } /// @brief Calculate the squared 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors. /// /// There is a scheme ambiguity for momentum-type four vectors /// as to whether the pseudorapidity (a purely geometric concept) or the /// rapidity (a relativistic energy-momentum quantity) is to be used: this can /// be chosen via the optional scheme parameter. inline double deltaR2(const FourVector& v, double eta2, double phi2, RapScheme scheme=PSEUDORAPIDITY) { switch (scheme) { case PSEUDORAPIDITY : return deltaR2(v.vector3(), eta2, phi2); case RAPIDITY: { const FourMomentum* mv = dynamic_cast(&v); if (!mv) { string err = "deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors"; throw std::runtime_error(err); } return deltaR2(*mv, eta2, phi2, scheme); } default: throw std::runtime_error("The specified deltaR scheme is not yet implemented"); } } /// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors. /// /// There is a scheme ambiguity for momentum-type four vectors /// as to whether the pseudorapidity (a purely geometric concept) or the /// rapidity (a relativistic energy-momentum quantity) is to be used: this can /// be chosen via the optional scheme parameter. inline double deltaR(const FourVector& v, double eta2, double phi2, RapScheme scheme=PSEUDORAPIDITY) { return sqrt(deltaR2(v, eta2, phi2, scheme)); } /// @brief Calculate the squared 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors. /// /// There is a scheme ambiguity for momentum-type four vectors /// as to whether the pseudorapidity (a purely geometric concept) or the /// rapidity (a relativistic energy-momentum quantity) is to be used: this can /// be chosen via the optional scheme parameter. inline double deltaR2(double eta1, double phi1, const FourVector& v, RapScheme scheme=PSEUDORAPIDITY) { switch (scheme) { case PSEUDORAPIDITY : return deltaR2(eta1, phi1, v.vector3()); case RAPIDITY: { const FourMomentum* mv = dynamic_cast(&v); if (!mv) { string err = "deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors"; throw std::runtime_error(err); } return deltaR2(eta1, phi1, *mv, scheme); } default: throw std::runtime_error("The specified deltaR scheme is not yet implemented"); } } /// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors. /// /// There is a scheme ambiguity for momentum-type four vectors /// as to whether the pseudorapidity (a purely geometric concept) or the /// rapidity (a relativistic energy-momentum quantity) is to be used: this can /// be chosen via the optional scheme parameter. inline double deltaR(double eta1, double phi1, const FourVector& v, RapScheme scheme=PSEUDORAPIDITY) { return sqrt(deltaR2(eta1, phi1, v, scheme)); } /// @brief Calculate the squared 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors. /// /// There is a scheme ambiguity for momentum-type four vectors /// as to whether the pseudorapidity (a purely geometric concept) or the /// rapidity (a relativistic energy-momentum quantity) is to be used: this can /// be chosen via the optional scheme parameter. inline double deltaR2(const FourMomentum& a, const FourMomentum& b, RapScheme scheme=PSEUDORAPIDITY) { switch (scheme) { case PSEUDORAPIDITY: return deltaR2(a.vector3(), b.vector3()); case RAPIDITY: return deltaR2(a.rapidity(), a.azimuthalAngle(), b.rapidity(), b.azimuthalAngle()); default: throw std::runtime_error("The specified deltaR scheme is not yet implemented"); } } /// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors. /// /// There is a scheme ambiguity for momentum-type four vectors /// as to whether the pseudorapidity (a purely geometric concept) or the /// rapidity (a relativistic energy-momentum quantity) is to be used: this can /// be chosen via the optional scheme parameter. inline double deltaR(const FourMomentum& a, const FourMomentum& b, RapScheme scheme=PSEUDORAPIDITY) { return sqrt(deltaR2(a, b, scheme)); } /// @brief Calculate the squared 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors. /// There is a scheme ambiguity for momentum-type four vectors /// as to whether the pseudorapidity (a purely geometric concept) or the /// rapidity (a relativistic energy-momentum quantity) is to be used: this can /// be chosen via the optional scheme parameter. inline double deltaR2(const FourMomentum& v, double eta2, double phi2, RapScheme scheme=PSEUDORAPIDITY) { switch (scheme) { case PSEUDORAPIDITY: return deltaR2(v.vector3(), eta2, phi2); case RAPIDITY: return deltaR2(v.rapidity(), v.azimuthalAngle(), eta2, phi2); default: throw std::runtime_error("The specified deltaR scheme is not yet implemented"); } } /// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors. /// There is a scheme ambiguity for momentum-type four vectors /// as to whether the pseudorapidity (a purely geometric concept) or the /// rapidity (a relativistic energy-momentum quantity) is to be used: this can /// be chosen via the optional scheme parameter. inline double deltaR(const FourMomentum& v, double eta2, double phi2, RapScheme scheme=PSEUDORAPIDITY) { return sqrt(deltaR2(v, eta2, phi2, scheme)); } /// @brief Calculate the squared 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors. /// There is a scheme ambiguity for momentum-type four vectors /// as to whether the pseudorapidity (a purely geometric concept) or the /// rapidity (a relativistic energy-momentum quantity) is to be used: this can /// be chosen via the optional scheme parameter. inline double deltaR2(double eta1, double phi1, const FourMomentum& v, RapScheme scheme=PSEUDORAPIDITY) { switch (scheme) { case PSEUDORAPIDITY: return deltaR2(eta1, phi1, v.vector3()); case RAPIDITY: return deltaR2(eta1, phi1, v.rapidity(), v.azimuthalAngle()); default: throw std::runtime_error("The specified deltaR scheme is not yet implemented"); } } /// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors. /// There is a scheme ambiguity for momentum-type four vectors /// as to whether the pseudorapidity (a purely geometric concept) or the /// rapidity (a relativistic energy-momentum quantity) is to be used: this can /// be chosen via the optional scheme parameter. inline double deltaR(double eta1, double phi1, const FourMomentum& v, RapScheme scheme=PSEUDORAPIDITY) { return sqrt(deltaR2(eta1, phi1, v, scheme)); } /// @brief Calculate the squared 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors. /// There is a scheme ambiguity for momentum-type four vectors /// as to whether the pseudorapidity (a purely geometric concept) or the /// rapidity (a relativistic energy-momentum quantity) is to be used: this can /// be chosen via the optional scheme parameter. inline double deltaR2(const FourMomentum& a, const FourVector& b, RapScheme scheme=PSEUDORAPIDITY) { switch (scheme) { case PSEUDORAPIDITY: return deltaR2(a.vector3(), b.vector3()); case RAPIDITY: return deltaR2(a.rapidity(), a.azimuthalAngle(), FourMomentum(b).rapidity(), b.azimuthalAngle()); default: throw std::runtime_error("The specified deltaR scheme is not yet implemented"); } } /// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors. /// There is a scheme ambiguity for momentum-type four vectors /// as to whether the pseudorapidity (a purely geometric concept) or the /// rapidity (a relativistic energy-momentum quantity) is to be used: this can /// be chosen via the optional scheme parameter. inline double deltaR(const FourMomentum& a, const FourVector& b, RapScheme scheme=PSEUDORAPIDITY) { return sqrt(deltaR2(a, b, scheme)); } /// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors. /// There is a scheme ambiguity for momentum-type four vectors /// as to whether the pseudorapidity (a purely geometric concept) or the /// rapidity (a relativistic energy-momentum quantity) is to be used: this can /// be chosen via the optional scheme parameter. inline double deltaR2(const FourVector& a, const FourMomentum& b, RapScheme scheme=PSEUDORAPIDITY) { return deltaR2(b, a, scheme); //< note reversed args } /// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors. /// There is a scheme ambiguity for momentum-type four vectors /// as to whether the pseudorapidity (a purely geometric concept) or the /// rapidity (a relativistic energy-momentum quantity) is to be used: this can /// be chosen via the optional scheme parameter. inline double deltaR(const FourVector& a, const FourMomentum& b, RapScheme scheme=PSEUDORAPIDITY) { return deltaR(b, a, scheme); //< note reversed args } /// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between a /// three-vector and a four-vector. inline double deltaR2(const FourMomentum& a, const Vector3& b) { return deltaR2(a.vector3(), b); } /// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between a /// three-vector and a four-vector. inline double deltaR(const FourMomentum& a, const Vector3& b) { return deltaR(a.vector3(), b); } /// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between a /// three-vector and a four-vector. inline double deltaR2(const Vector3& a, const FourMomentum& b) { return deltaR2(a, b.vector3()); } /// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between a /// three-vector and a four-vector. inline double deltaR(const Vector3& a, const FourMomentum& b) { return deltaR(a, b.vector3()); } /// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between a /// three-vector and a four-vector. inline double deltaR2(const FourVector& a, const Vector3& b) { return deltaR2(a.vector3(), b); } /// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between a /// three-vector and a four-vector. inline double deltaR(const FourVector& a, const Vector3& b) { return deltaR(a.vector3(), b); } /// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between a /// three-vector and a four-vector. inline double deltaR2(const Vector3& a, const FourVector& b) { return deltaR2(a, b.vector3()); } /// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between a /// three-vector and a four-vector. inline double deltaR(const Vector3& a, const FourVector& b) { return deltaR(a, b.vector3()); } //@} ////////////////////////////////////////////////////// /// @name \f$ \Delta phi \f$ calculations from 4-vectors //@{ /// Calculate the difference in azimuthal angle between two vectors. inline double deltaPhi(const FourMomentum& a, const FourMomentum& b) { return deltaPhi(a.vector3(), b.vector3()); } /// Calculate the difference in azimuthal angle between two vectors. inline double deltaPhi(const FourMomentum& v, double phi2) { return deltaPhi(v.vector3(), phi2); } /// Calculate the difference in azimuthal angle between two vectors. inline double deltaPhi(double phi1, const FourMomentum& v) { return deltaPhi(phi1, v.vector3()); } /// Calculate the difference in azimuthal angle between two vectors. inline double deltaPhi(const FourVector& a, const FourVector& b) { return deltaPhi(a.vector3(), b.vector3()); } /// Calculate the difference in azimuthal angle between two vectors. inline double deltaPhi(const FourVector& v, double phi2) { return deltaPhi(v.vector3(), phi2); } /// Calculate the difference in azimuthal angle between two vectors. inline double deltaPhi(double phi1, const FourVector& v) { return deltaPhi(phi1, v.vector3()); } /// Calculate the difference in azimuthal angle between two vectors. inline double deltaPhi(const FourVector& a, const FourMomentum& b) { return deltaPhi(a.vector3(), b.vector3()); } /// Calculate the difference in azimuthal angle between two vectors. inline double deltaPhi(const FourMomentum& a, const FourVector& b) { return deltaPhi(a.vector3(), b.vector3()); } /// Calculate the difference in azimuthal angle between two vectors. inline double deltaPhi(const FourVector& a, const Vector3& b) { return deltaPhi(a.vector3(), b); } /// Calculate the difference in azimuthal angle between two vectors. inline double deltaPhi(const Vector3& a, const FourVector& b) { return deltaPhi(a, b.vector3()); } /// Calculate the difference in azimuthal angle between two vectors. inline double deltaPhi(const FourMomentum& a, const Vector3& b) { return deltaPhi(a.vector3(), b); } /// Calculate the difference in azimuthal angle between two vectors. inline double deltaPhi(const Vector3& a, const FourMomentum& b) { return deltaPhi(a, b.vector3()); } //@} ////////////////////////////////////////////////////// /// @name \f$ |\Delta eta| \f$ calculations from 4-vectors //@{ /// Calculate the difference in pseudorapidity between two vectors. inline double deltaEta(const FourMomentum& a, const FourMomentum& b) { return deltaEta(a.vector3(), b.vector3()); } /// Calculate the difference in pseudorapidity between two vectors. inline double deltaEta(const FourMomentum& v, double eta2) { return deltaEta(v.vector3(), eta2); } /// Calculate the difference in pseudorapidity between two vectors. inline double deltaEta(double eta1, const FourMomentum& v) { return deltaEta(eta1, v.vector3()); } /// Calculate the difference in pseudorapidity between two vectors. inline double deltaEta(const FourVector& a, const FourVector& b) { return deltaEta(a.vector3(), b.vector3()); } /// Calculate the difference in pseudorapidity between two vectors. inline double deltaEta(const FourVector& v, double eta2) { return deltaEta(v.vector3(), eta2); } /// Calculate the difference in pseudorapidity between two vectors. inline double deltaEta(double eta1, const FourVector& v) { return deltaEta(eta1, v.vector3()); } /// Calculate the difference in pseudorapidity between two vectors. inline double deltaEta(const FourVector& a, const FourMomentum& b) { return deltaEta(a.vector3(), b.vector3()); } /// Calculate the difference in pseudorapidity between two vectors. inline double deltaEta(const FourMomentum& a, const FourVector& b) { return deltaEta(a.vector3(), b.vector3()); } /// Calculate the difference in pseudorapidity between two vectors. inline double deltaEta(const FourVector& a, const Vector3& b) { return deltaEta(a.vector3(), b); } /// Calculate the difference in pseudorapidity between two vectors. inline double deltaEta(const Vector3& a, const FourVector& b) { return deltaEta(a, b.vector3()); } /// Calculate the difference in pseudorapidity between two vectors. inline double deltaEta(const FourMomentum& a, const Vector3& b) { return deltaEta(a.vector3(), b); } /// Calculate the difference in pseudorapidity between two vectors. inline double deltaEta(const Vector3& a, const FourMomentum& b) { return deltaEta(a, b.vector3()); } //@} /// @name \f$ |\Delta y| \f$ calculations from 4-momentum vectors //@{ /// Calculate the difference in rapidity between two 4-momentum vectors. inline double deltaRap(const FourMomentum& a, const FourMomentum& b) { return deltaRap(a.rapidity(), b.rapidity()); } /// Calculate the difference in rapidity between two 4-momentum vectors. inline double deltaRap(const FourMomentum& v, double y2) { return deltaRap(v.rapidity(), y2); } /// Calculate the difference in rapidity between two 4-momentum vectors. inline double deltaRap(double y1, const FourMomentum& v) { return deltaRap(y1, v.rapidity()); } //@} ////////////////////////////////////////////////////// /// @name 4-vector comparison functions (for sorting) //@{ /// Comparison to give a sorting by decreasing pT inline bool cmpMomByPt(const FourMomentum& a, const FourMomentum& b) { return a.pt() > b.pt(); } /// Comparison to give a sorting by increasing pT inline bool cmpMomByAscPt(const FourMomentum& a, const FourMomentum& b) { return a.pt() < b.pt(); } /// Comparison to give a sorting by decreasing 3-momentum magnitude |p| inline bool cmpMomByP(const FourMomentum& a, const FourMomentum& b) { return a.vector3().mod() > b.vector3().mod(); } /// Comparison to give a sorting by increasing 3-momentum magnitude |p| inline bool cmpMomByAscP(const FourMomentum& a, const FourMomentum& b) { return a.vector3().mod() < b.vector3().mod(); } /// Comparison to give a sorting by decreasing transverse energy inline bool cmpMomByEt(const FourMomentum& a, const FourMomentum& b) { return a.Et() > b.Et(); } /// Comparison to give a sorting by increasing transverse energy inline bool cmpMomByAscEt(const FourMomentum& a, const FourMomentum& b) { return a.Et() < b.Et(); } /// Comparison to give a sorting by decreasing energy inline bool cmpMomByE(const FourMomentum& a, const FourMomentum& b) { return a.E() > b.E(); } /// Comparison to give a sorting by increasing energy inline bool cmpMomByAscE(const FourMomentum& a, const FourMomentum& b) { return a.E() < b.E(); } /// Comparison to give a sorting by decreasing mass inline bool cmpMomByMass(const FourMomentum& a, const FourMomentum& b) { return a.mass() > b.mass(); } /// Comparison to give a sorting by increasing mass inline bool cmpMomByAscMass(const FourMomentum& a, const FourMomentum& b) { return a.mass() < b.mass(); } /// Comparison to give a sorting by increasing eta (pseudorapidity) inline bool cmpMomByEta(const FourMomentum& a, const FourMomentum& b) { return a.eta() < b.eta(); } /// Comparison to give a sorting by decreasing eta (pseudorapidity) inline bool cmpMomByDescEta(const FourMomentum& a, const FourMomentum& b) { return a.pseudorapidity() > b.pseudorapidity(); } /// Comparison to give a sorting by increasing absolute eta (pseudorapidity) inline bool cmpMomByAbsEta(const FourMomentum& a, const FourMomentum& b) { return fabs(a.eta()) < fabs(b.eta()); } /// Comparison to give a sorting by increasing absolute eta (pseudorapidity) inline bool cmpMomByDescAbsEta(const FourMomentum& a, const FourMomentum& b) { return fabs(a.eta()) > fabs(b.eta()); } /// Comparison to give a sorting by increasing rapidity inline bool cmpMomByRap(const FourMomentum& a, const FourMomentum& b) { return a.rapidity() < b.rapidity(); } /// Comparison to give a sorting by decreasing rapidity inline bool cmpMomByDescRap(const FourMomentum& a, const FourMomentum& b) { return a.rapidity() > b.rapidity(); } /// Comparison to give a sorting by increasing absolute rapidity inline bool cmpMomByAbsRap(const FourMomentum& a, const FourMomentum& b) { return fabs(a.rapidity()) < fabs(b.rapidity()); } /// Comparison to give a sorting by decreasing absolute rapidity inline bool cmpMomByDescAbsRap(const FourMomentum& a, const FourMomentum& b) { return fabs(a.rapidity()) > fabs(b.rapidity()); } /// @todo Add sorting by phi [0..2PI] /// Sort a container of momenta by cmp and return by reference for non-const inputs template inline MOMS& sortBy(MOMS& pbs, const CMP& cmp) { std::sort(pbs.begin(), pbs.end(), cmp); return pbs; } /// Sort a container of momenta by cmp and return by value for const inputs template inline MOMS sortBy(const MOMS& pbs, const CMP& cmp) { MOMS rtn = pbs; std::sort(rtn.begin(), rtn.end(), cmp); return rtn; } /// Sort a container of momenta by pT (decreasing) and return by reference for non-const inputs template inline MOMS& sortByPt(MOMS& pbs) { return sortBy(pbs, cmpMomByPt); } /// Sort a container of momenta by pT (decreasing) and return by value for const inputs template inline MOMS sortByPt(const MOMS& pbs) { return sortBy(pbs, cmpMomByPt); } /// Sort a container of momenta by E (decreasing) and return by reference for non-const inputs template inline MOMS& sortByE(MOMS& pbs) { return sortBy(pbs, cmpMomByE); } /// Sort a container of momenta by E (decreasing) and return by value for const inputs template inline MOMS sortByE(const MOMS& pbs) { return sortBy(pbs, cmpMomByE); } /// Sort a container of momenta by Et (decreasing) and return by reference for non-const inputs template inline MOMS& sortByEt(MOMS& pbs) { return sortBy(pbs, cmpMomByEt); } /// Sort a container of momenta by Et (decreasing) and return by value for const inputs template inline MOMS sortByEt(const MOMS& pbs) { return sortBy(pbs, cmpMomByEt); } //@} /// @name MT calculation //@{ /// Calculate transverse mass of a visible and an invisible 3-vector inline double mT(const Vector3& vis, const Vector3& invis) { // return sqrt(2*vis.perp()*invis.perp() * (1 - cos(deltaPhi(vis, invis))) ); return mT(vis.perp(), invis.perp(), deltaPhi(vis, invis)); } /// Calculate transverse mass of a visible and an invisible 4-vector inline double mT(const FourMomentum& vis, const FourMomentum& invis) { return mT(vis.p3(), invis.p3()); } /// Calculate transverse mass of a visible 4-vector and an invisible 3-vector inline double mT(const FourMomentum& vis, const Vector3& invis) { return mT(vis.p3(), invis); } /// Calculate transverse mass of a visible 4-vector and an invisible 3-vector inline double mT(const Vector3& vis, const FourMomentum& invis) { return mT(vis, invis.p3()); } //@} ////////////////////////////////////////////////////// /// @name 4-vector string representations //@{ /// Render a 4-vector as a string. inline std::string toString(const FourVector& lv) { std::ostringstream out; out << "(" << (fabs(lv.t()) < 1E-30 ? 0.0 : lv.t()) << "; " << (fabs(lv.x()) < 1E-30 ? 0.0 : lv.x()) << ", " << (fabs(lv.y()) < 1E-30 ? 0.0 : lv.y()) << ", " << (fabs(lv.z()) < 1E-30 ? 0.0 : lv.z()) << ")"; return out.str(); } /// Write a 4-vector to an ostream. inline std::ostream& operator<<(std::ostream& out, const FourVector& lv) { out << toString(lv); return out; } //@} /// @name Typedefs for lists of vector types //@{ typedef std::vector FourVectors; typedef std::vector FourMomenta; //@} } #endif diff --git a/src/Projections/DISKinematics.cc b/src/Projections/DISKinematics.cc --- a/src/Projections/DISKinematics.cc +++ b/src/Projections/DISKinematics.cc @@ -1,80 +1,80 @@ // -*- C++ -*- #include "Rivet/Projections/DISKinematics.hh" #include "Rivet/Math/Constants.hh" namespace Rivet { void DISKinematics::project(const Event& e) { // Find appropriate DIS leptons const DISLepton& dislep = applyProjection(e, "Lepton"); _outLepton = dislep.out(); // Identify beam hadron const ParticlePair& inc = applyProjection(e, "Beam").beams(); const bool firstIsHadron = PID::isHadron(inc.first.pid()); const bool secondIsHadron = PID::isHadron(inc.second.pid()); if (firstIsHadron && !secondIsHadron) { _inHadron = inc.first; _inLepton = dislep.in(); // inc.second; } else if (!firstIsHadron && secondIsHadron) { _inHadron = inc.second; _inLepton = dislep.in(); // inc.first; } else { throw Error("DISKinematics could not find the correct beam hadron"); } // Get the DIS lepton and store some of its properties const FourMomentum pHad = _inHadron.momentum(); const FourMomentum pLepIn = _inLepton.momentum(); const FourMomentum pLepOut = _outLepton.momentum(); const FourMomentum pGamma = pLepIn - pLepOut; const FourMomentum tothad = pGamma + pHad; _theQ2 = -pGamma.mass2(); _theW2 = tothad.mass2(); _theX = Q2()/(2.0 * pGamma * pHad); _theY = (pGamma * pHad) / (pLepIn * pHad); _theS = invariant(pLepIn + pHad); // Calculate boost vector for boost into HCM-system LorentzTransform tmp; - tmp.setBetaVec(-tothad.boostVector()); + tmp.setBetaVec(-tothad.betaVec()); // Rotate so the photon is in x-z plane in HCM rest frame FourMomentum pGammaHCM = tmp.transform(pGamma); tmp.preMult(Matrix3(Vector3::mkZ(), -pGammaHCM.azimuthalAngle())); pGammaHCM = tmp.transform(pGamma); assert(isZero(dot(pGammaHCM.vector3(), Vector3::mkY()))); // Rotate so the photon is along the positive z-axis const double rot_angle = pGammaHCM.polarAngle() * (pGammaHCM.px() >= 0 ? -1 : 1); tmp.preMult(Matrix3(Vector3::mkY(), rot_angle)); // Check that final HCM photon lies along +ve z as expected pGammaHCM = tmp.transform(pGamma); assert(isZero(dot(pGammaHCM.vector3(), Vector3::mkX()), 1e-3)); assert(isZero(dot(pGammaHCM.vector3(), Vector3::mkY()), 1e-3)); assert(isZero(angle(pGammaHCM.vector3(), Vector3::mkZ()), 1e-3)); // Finally rotate so outgoing lepton at phi = 0 FourMomentum pLepOutHCM = tmp.transform(pLepOut); tmp.preMult(Matrix3(Vector3::mkZ(), -pLepOutHCM.azimuthalAngle())); assert(isZero(tmp.transform(pLepOut).azimuthalAngle())); _hcm = tmp; // Boost to Breit frame (use opposite convention for photon --- along *minus* z) tmp.preMult(Matrix3(Vector3::mkX(), PI)); const double bz = 1 - 2*x(); _breit = LorentzTransform::mkObjTransformFromBeta(Vector3::mkZ() * bz).combine(tmp); assert(isZero(angle(_breit.transform(pGamma).vector3(), -Vector3::mkZ()), 1e-3)); assert(isZero(_breit.transform(pLepOut).azimuthalAngle(), 1e-3)); } CmpState DISKinematics::compare(const Projection & p) const { const DISKinematics& other = pcast(p); return mkNamedPCmp(other, "Lepton"); } } diff --git a/test/testMatVec.cc b/test/testMatVec.cc --- a/test/testMatVec.cc +++ b/test/testMatVec.cc @@ -1,155 +1,155 @@ #include "Rivet/Math/MathUtils.hh" #include "Rivet/Math/Vectors.hh" #include "Rivet/Math/Matrices.hh" // #include "Rivet/Math/MatrixDiag.hh" using namespace Rivet; #include #include #include using namespace std; int main() { FourVector a(1,0,0,0); cout << a << ": interval = " << a.invariant() << '\n'; assert(fuzzyEquals(a.invariant(), 1)); a.setZ(1); assert(isZero(a.invariant())); cout << a << ": interval = " << a.invariant() << '\n'; a.setY(2).setZ(3); cout << a << ": interval = " << a.invariant() << '\n'; assert(fuzzyEquals(a.invariant(), -12)); cout << a << ": vector = " << a.vector3() << 'n' << '\n'; FourMomentum b(1,0,0,0); cout << b << ": mass = " << b.mass() << '\n'; assert(fuzzyEquals(b.mass2(), 1)); b.setPz(1); cout << b << ": mass = " << b.mass() << '\n'; assert(isZero(b.mass2())); b.setPy(2).setPz(3).setE(6); cout << b << ": mass = " << b.mass2() << '\n'; assert(fuzzyEquals(b.mass2(), 23)); cout << b << ": vector = " << b.vector3() << 'n' << '\n'; Matrix3 m; m.set(0, 0, 7/4.0); m.set(0, 1, 3 * sqrt(3)/4.0); m.set(1, 0, 3 * sqrt(3)/4.0); m.set(1, 1, 13/4.0); m.set(2, 2, 9); cout << m << 'n' << '\n'; // EigenSystem<3> es = diagonalize(m); cout << "Matrices:" << '\n'; cout << Matrix3() << '\n'; cout << Matrix3::mkIdentity() << '\n'; const Matrix3 I3 = Matrix3::mkIdentity(); cout << Matrix3::mkIdentity() * m * I3 << '\n'; cout << "tr(0) & det(0): " << Matrix3().trace() << ", " << Matrix3().det() << '\n'; cout << "tr(I3) & det(I3): " << I3.trace() << ", " << I3.det() << '\n'; Matrix3 m1 = Matrix3::mkIdentity(); Matrix3 m2 = m1; m1.setRow(1, Vector3(1,2,3)); m2.setColumn(1, Vector3(3,2,1)); Matrix3 m3 = Matrix3::mkZero(); cout << m1 << " + " << m2 << " = " << m1 + m2 << '\n'; m3.setRow(0, Vector3(2,3,0)).setRow(1, Vector3(1,4,3)).setRow(2, Vector3(0,1,2)); cout << m1+m2 << " == " << m3 << ": " << (m1+m2 == m3 ? "true" : "false") << '\n'; cout << '\n'; Vector3 v3(1,2,3); cout << "Vector: " << v3 << '\n'; cout << "Invert: " << v3 << " --> " << -v3 << '\n'; const Matrix3 rot90(Vector3(0,0,1), PI/2.0); const Matrix3 rot90m(Vector3(0,0,1), -PI/2.0); const Matrix3 rot180(Vector3(0,0,1), PI); const Matrix3 rot180m(Vector3(0,0,1), -PI); const Vector3 v3_90 = rot90*v3; cout << "Rot 90: " << v3 << " ---> " << v3_90 << '\n'; const Vector3 v3_90m = rot90m*v3; cout << "Rot -90: " << v3 << " ---> " << v3_90m << '\n'; const Vector3 v3_180 = rot180*v3; cout << "Rot 180: " << v3 << " ---> " << v3_180 << '\n'; const Vector3 v3_180m = rot180m*v3; cout << "Rot -180: " << v3 << " ---> " << v3_180m << '\n'; assert(fuzzyEquals(v3_180, v3_180m)); const Vector3 v3_9090 = rot90*rot90*v3; cout << "Rot 2 x 90: " << v3 << " ---> " << v3_9090 << '\n'; assert(fuzzyEquals(v3_180, v3_9090)); const Vector3 v3_90m90m = rot90m*rot90m*v3; cout << "Rot 2 x -90: " << v3 << " ---> " << v3_90m90m << '\n'; assert(fuzzyEquals(v3_180, v3_90m90m)); const Vector3 v3_9090m = rot90*rot90m*v3; const Vector3 v3_90m90 = rot90m*rot90*v3; cout << "Rot 90*-90: "<< v3 << " ---> " << v3_9090m << '\n'; cout << "Rot -90*90: "<< v3 << " ---> " << v3_90m90 << '\n'; assert(fuzzyEquals(v3, v3_9090m)); assert(fuzzyEquals(v3, v3_90m90)); const Vector3 v3_90i = rot90.inverse()*v3; cout << "Rot (90)^-1: "<< v3 << " ---> " << v3_90i << '\n'; assert(fuzzyEquals(v3_90i, v3_90m)); const Vector3 v3_9090i = rot90*rot90.inverse()*v3; const Vector3 v3_90i90 = rot90.inverse()*rot90*v3; cout << "Rot 90*(90)^-1: "<< v3 << " ---> " << v3_9090i << '\n'; cout << "Rot (90)^-1*90: "<< v3 << " ---> " << v3_90i90 << '\n'; assert(fuzzyEquals(v3, v3_9090i)); assert(fuzzyEquals(v3, v3_90i90)); const Matrix3 rot1(Vector3(0,1,0), PI/180.0); cout << "Rot 0 x 45 x 1: " << v3 << '\n'; for (size_t i = 0; i < 8; ++i) { for (size_t j = 0; j < 45; ++j) { v3 = rot1*v3; } cout << "Rot " << i+1 << " x 45 x 1: " << v3 << '\n'; } assert(fuzzyEquals(v3, Vector3(1,2,3))); cout << '\n'; cout << "Boosts:" << '\n'; LorentzTransform ltX = LorentzTransform::mkObjTransformFromBeta(Vector3(0.5, 0, 0)); cout << "LTx: " << ltX << '\n'; cout << "I on LTx: " << ltX.rotate(Matrix3::mkIdentity()) << '\n'; cout << "Rot90 on LTx: " << ltX.rotate(rot90) << '\n'; cout << '\n'; cout << "X-boosts:" << '\n'; const FourMomentum p1 = FourMomentum(10, 0, 0, 1); const FourMomentum p2 = ltX.transform(p1); cout << p1 << " -> " << p2 << '\n'; cout << p2 << " -> " << ltX.inverse().transform(p2) << '\n'; - //cout << p1.boostVector() << '\n'; - const FourMomentum p3 = LorentzTransform::mkFrameTransformFromBeta(p1.boostVector()).transform(p1); + //cout << p1.betaVec() << '\n'; + const FourMomentum p3 = LorentzTransform::mkFrameTransformFromBeta(p1.betaVec()).transform(p1); cout << p1 << " -> " << p3 << '\n'; cout << '\n'; LorentzTransform ltY = LorentzTransform::mkObjTransformFromGamma(Vector3(0, 1.2, 0)); cout << FourMomentum(1,0,0,1) << " -> " //<< "\n " << (ltX * ltY).transform(FourMomentum(1,0,0,1)) << '\n'; cout << FourMomentum(1,0,0,1) << " -> " //<< "\n " << (ltY * ltX).transform(FourMomentum(1,0,0,1)) << '\n'; cout << (ltX * ltY).betaVec() << '\n'; cout << (ltY * ltX).betaVec() << '\n'; cout << (ltX * ltX.inverse()).betaVec() << '\n'; // If we are already in the rest frame and there is no boost, then LT is trivial/identity LorentzTransform noBoost; cout << "Element 0,0 should be 1 and is " << noBoost.toMatrix().get(0,0) << '\n'; assert(noBoost.toMatrix().get(0,0)==1); cout << "Element 0,1 should be 0 and is " << noBoost.toMatrix().get(0,1) << '\n'; assert(noBoost.toMatrix().get(0,1)==0); cout << "Element 1,0 should be 0 and is " << noBoost.toMatrix().get(1,0) << '\n'; assert(noBoost.toMatrix().get(1,0)==0); cout << "Element 1,1 should be 1 and is " << noBoost.toMatrix().get(1,1) << '\n'; assert(noBoost.toMatrix().get(1,1)==1); return EXIT_SUCCESS; }