diff --git a/analyses/pluginALICE/ALICE_2012_I930312.cc b/analyses/pluginALICE/ALICE_2012_I930312.cc --- a/analyses/pluginALICE/ALICE_2012_I930312.cc +++ b/analyses/pluginALICE/ALICE_2012_I930312.cc @@ -1,335 +1,316 @@ // -*- C++ -*- #include "Rivet/Projections/ChargedFinalState.hh" -#include "Rivet/Tools/Cuts.hh" #include "Rivet/Projections/SingleValueProjection.hh" #include "Rivet/Tools/AliceCommon.hh" #include "Rivet/Projections/AliceCommon.hh" -#include - -#define _USE_MATH_DEFINES -#include namespace Rivet { - - /// Analysis + + class ALICE_2012_I930312 : public Analysis { - public: - + + // Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ALICE_2012_I930312); - + + /// Book histograms and initialise projections before the run void init() { // Declare centrality projection - declareCentrality(ALICE::V0MMultiplicity(), "ALICE_2015_PBPBCentrality", - "V0M", "V0M"); + declareCentrality(ALICE::V0MMultiplicity(), "ALICE_2015_PBPBCentrality", "V0M", "V0M"); - // Charged final states with |eta| < 1.0 and 8 < pT < 15 GeV/c for - // trigger particles - const Cut& cutTrigger = - Cuts::abseta < 1.0 && Cuts::pT > 8*GeV && Cuts::pT < 15*GeV; + // Charged final states with |eta| < 1.0 and 8 < pT < 15 GeV/c for trigger particles + const Cut& cutTrigger = Cuts::abseta < 1.0 && Cuts::pT > 8*GeV && Cuts::pT < 15*GeV; const ChargedFinalState cfsTrigger(cutTrigger); addProjection(cfsTrigger,"CFSTrigger"); - + // Set limit values of pT bins pt_limits[0] = 3.; pt_limits[1] = 4.; pt_limits[2] = 6.; pt_limits[3] = 8.; pt_limits[4] = 10.; - - // Charged final states with |eta| < 1.0 and different pT bins for - // associated particles. + + // Charged final states with |eta| < 1.0 and different pT bins for associated particles. for (int ipt = 0; ipt < PT_BINS; ipt++) { - Cut mycut = Cuts::abseta < 1.0 && Cuts::pT > pt_limits[ipt]*GeV && - Cuts::pT < pt_limits[ipt + 1]*GeV; - declare(ChargedFinalState(mycut), "CFSAssoc" + std::to_string(ipt)); + Cut mycut = Cuts::abseta < 1.0 && Cuts::pT > pt_limits[ipt]*GeV && Cuts::pT < pt_limits[ipt + 1]*GeV; + declare(ChargedFinalState(mycut), "CFSAssoc" + std::to_string(ipt)); } - + // Create event strings event_string[0] = "pp"; event_string[1] = "central"; event_string[2] = "peripheral"; event_string[3] = "other"; - + // For each event type for (int itype = 0; itype < EVENT_TYPES; itype++) { - // For each pT range - for (int ipt = 0; ipt < PT_BINS; ipt++) { - // Initialize yield histograms - _histYield[itype][ipt] = bookHisto1D("Yield_" + event_string[itype] - + "_" + std::to_string(ipt), 36, -0.5 * M_PI, 1.5 * M_PI, - "Associated particle per trigger particle yield", - "$\\Delta\\eta$ (rad)", "$1 / N_{trig} dN_{assoc} / d\\Delta\\eta$ (rad$^-1$)"); - _histYieldBkgRemoved[itype][ipt] = bookHisto1D("Yield_" + - event_string[itype] + "_nobkg_" + std::to_string(ipt), 36, -0.5*M_PI, - 1.5 * M_PI, "Associated particle per trigger particle yield no bkg", - "$\\Delta\\eta$ (rad)", "$1 / N_{trig} dN_{assoc} / d\\Delta\\eta$ (rad$^-1$)"); - } + // For each pT range + for (int ipt = 0; ipt < PT_BINS; ipt++) { + // Initialize yield histograms + _histYield[itype][ipt] = bookHisto1D("Yield_" + event_string[itype] + + "_" + std::to_string(ipt), 36, -0.5 * M_PI, 1.5 * M_PI, + "Associated particle per trigger particle yield", + "$\\Delta\\eta$ (rad)", "$1 / N_{trig} dN_{assoc} / d\\Delta\\eta$ (rad$^-1$)"); + _histYieldBkgRemoved[itype][ipt] = bookHisto1D("Yield_" + + event_string[itype] + "_nobkg_" + std::to_string(ipt), 36, -0.5*M_PI, + 1.5 * M_PI, "Associated particle per trigger particle yield no bkg", + "$\\Delta\\eta$ (rad)", "$1 / N_{trig} dN_{assoc} / d\\Delta\\eta$ (rad$^-1$)"); + } } - + // Histogram for counting trigger particles for each event type - _histTriggerCounter = bookHisto1D("Trigger", EVENT_TYPES, 0.0, - EVENT_TYPES, "Trigger counter", "event type", "N"); - + _histTriggerCounter = bookHisto1D("Trigger", EVENT_TYPES, 0.0, + EVENT_TYPES, "Trigger counter", "event type", "N"); + } + /// Perform the per-event analysis void analyze(const Event& event) { - + const double weight = event.weight(); - + // Create charged final state for trigger particle - const ChargedFinalState& triggerFinalState = - applyProjection(event, "CFSTrigger"); + const ChargedFinalState& triggerFinalState = apply(event, "CFSTrigger"); Particles triggerParticles = triggerFinalState.particlesByPt(); - // Create charged final state for associated particle + // Create charged final state for associated particle ChargedFinalState associatedFinalState[PT_BINS]; Particles associatedParticles[PT_BINS]; for (int ipt = 0; ipt < PT_BINS; ipt++) { - associatedFinalState[ipt] = applyProjection(event, - "CFSAssoc" + std::to_string(ipt)); - associatedParticles[ipt] = associatedFinalState[ipt].particlesByPt(); + associatedFinalState[ipt] = apply(event, "CFSAssoc" + std::to_string(ipt)); + associatedParticles[ipt] = associatedFinalState[ipt].particlesByPt(); } - + // Check event type if (event.genEvent()->heavy_ion()) { - // Prepare centrality projection and value - const CentralityProjection& centrProj = - apply(event, "V0M"); - double centr = centrProj(); - - // Set the flag for the type of the event - if (centr > 0.0 && centr < 5.0) - event_type = 1; // PbPb, central - else if (centr > 60.0 && centr < 90.0) - event_type = 2; // PbPb, peripherial - else - event_type = 3; // PbPb, other + // Prepare centrality projection and value + const CentralityProjection& centrProj = apply(event, "V0M"); + double centr = centrProj(); + + // Set the flag for the type of the event + if (centr > 0.0 && centr < 5.0) + event_type = 1; // PbPb, central + else if (centr > 60.0 && centr < 90.0) + event_type = 2; // PbPb, peripherial + else + event_type = 3; // PbPb, other } else { - event_type = 0; // pp + event_type = 0; // pp } - + // Veto event if not valid event type - if (event_type == 3) - vetoEvent; - + if (event_type == 3) vetoEvent; + // Fill trigger histogram for a proper event type _histTriggerCounter->fill(event_type, triggerParticles.size()); - + // Loop over trigger particles - for (const auto& triggerParticle : triggerParticles) { - // For each pt bin - for (int ipt = 0; ipt < PT_BINS; ipt++) { - // Loop over associated particles - for (const auto& associatedParticle : associatedParticles[ipt]) { - // If associated and trigger particle are not the same particles. - if (associatedParticle != triggerParticle) { - // Test trigger particle. - if (triggerParticle.pt() > associatedParticle.pt()) { - // Calculate delta phi in range (-0.5*PI, 1.5*PI). - double dPhi = triggerParticle.phi() - - associatedParticle.phi(); - while (dPhi > 1.5 * M_PI) { dPhi -= 2 * M_PI; } - while (dPhi < -0.5 * M_PI) { dPhi += 2 * M_PI; } - // Fill yield histogram for calculated delta phi - _histYield[event_type][ipt]->fill(dPhi, weight); - } - } - } - } + for (const Particle& triggerParticle : triggerParticles) { + // For each pt bin + for (int ipt = 0; ipt < PT_BINS; ipt++) { + // Loop over associated particles + for (const Particle& associatedParticle : associatedParticles[ipt]) { + // If associated and trigger particle are not the same particles. + if (!isSame(triggerParticle, associatedParticle)) { + // Test trigger particle. + if (triggerParticle.pt() > associatedParticle.pt()) { + // Calculate delta phi in range (-0.5*PI, 1.5*PI). + const double dPhi = deltaPhi(triggerParticle, associatedParticle, true); + // Fill yield histogram for calculated delta phi + _histYield[event_type][ipt]->fill(dPhi, weight); + } + } + } + } } } + /// Normalise histograms etc., after the run void finalize() { - + // Check for the reentrant finalize bool pp_available = false, PbPb_available = false; reentrant_flag = false; // For each event type for (int itype = 0; itype < EVENT_TYPES; itype++) { - // For each pT range - for (int ipt = 0; ipt < PT_BINS; ipt++) { - if (_histYield[itype][ipt]->numEntries() > 0) - itype == 0 ? pp_available = true : PbPb_available = true; - } + // For each pT range + for (int ipt = 0; ipt < PT_BINS; ipt++) { + if (_histYield[itype][ipt]->numEntries() > 0) + itype == 0 ? pp_available = true : PbPb_available = true; + } } reentrant_flag = (pp_available && PbPb_available); - + // Postprocessing of the histograms if (reentrant_flag == true) { - + // Initialize IAA and ICP histograms _histIAA[0] = bookScatter2D(1, 1, 1); _histIAA[1] = bookScatter2D(2, 1, 1); _histIAA[2] = bookScatter2D(5, 1, 1); - _histIAA[3] = bookScatter2D(3, 1, 1); _histIAA[4] = bookScatter2D(4, 1, 1); _histIAA[5] = bookScatter2D(6, 1, 1); - // Variables for near and away side peak calculation - double nearSide[EVENT_TYPES][PT_BINS] = { {0.0} }; - double awaySide[EVENT_TYPES][PT_BINS] = { {0.0} }; - - // Variables for background error calculation - double background[EVENT_TYPES][PT_BINS] = { {0.0} }; - double backgroundError[EVENT_TYPES][PT_BINS] = { {0.0} }; - - // Variables for integration error calculation - double scalingFactor[EVENT_TYPES] = {0.0}; - double numberOfEntries[EVENT_TYPES][PT_BINS][2] = { { {0.0} } }; - int numberOfBins[EVENT_TYPES][PT_BINS][2] = { { {0} } }; - - // For each event type - for (int itype = 0; itype < EVENT_TYPES; itype++) { - // For each pT range - for (int ipt = 0; ipt < PT_BINS; ipt++) { - // Check if histograms are fine - if (_histTriggerCounter->numEntries() == 0 || - _histYield[itype][ipt]->numEntries() == 0) { - cout << "There are no entries in one of the histograms" << endl; - continue; - } - - // Scale yield histogram - if ((_histTriggerCounter->bin(itype).sumW() != 0)) { - scalingFactor[itype] = 1. / - _histTriggerCounter->bin(itype).sumW(); - scale(_histYield[itype][ipt], - (1. / _histTriggerCounter->bin(itype).sumW())); - } - - // Calculate background - double sum = 0.0; - int nbins = 0; - for (unsigned int ibin = 0; ibin < _histYield[itype][ipt]->numBins(); ibin++) { - if ((_histYield[itype][ipt]->bin(ibin).xMid() > (-0.5 * M_PI) && - _histYield[itype][ipt]->bin(ibin).xMid() < (-0.5 * M_PI + 0.4)) || - (_histYield[itype][ipt]->bin(ibin).xMid() > (0.5 * M_PI - 0.4) && - _histYield[itype][ipt]->bin(ibin).xMid() < (0.5 * M_PI + 0.4)) || - (_histYield[itype][ipt]->bin(ibin).xMid() > (1.5 * M_PI - 0.4) && - _histYield[itype][ipt]->bin(ibin).xMid() < (1.5 * M_PI))) { - sum += _histYield[itype][ipt]->bin(ibin).sumW(); - nbins++; - } - } - if (nbins == 0) { - std::cout << "Failed to estimate background!" << std::endl; - continue; - } - background[itype][ipt] = sum / nbins; - - // Calculate background error - sum = 0.0; - nbins = 0; - for (unsigned int ibin = 0; ibin < _histYield[itype][ipt]->numBins(); ibin++) { - if (_histYield[itype][ipt]->bin(ibin).xMid() > (0.5 * M_PI - 0.4) && - _histYield[itype][ipt]->bin(ibin).xMid() < (0.5 * M_PI + 0.4)) { - sum += (_histYield[itype][ipt]->bin(ibin).sumW() - background[itype][ipt]) * - (_histYield[itype][ipt]->bin(ibin).sumW() - background[itype][ipt]); - nbins++; - } - } - backgroundError[itype][ipt] = sqrt(sum / (nbins - 1)); - - // Fill histograms with removed background - for (unsigned int ibin = 0; ibin < _histYield[itype][ipt]->numBins(); ibin++) { - _histYieldBkgRemoved[itype][ipt]->fillBin(ibin, - _histYield[itype][ipt]->bin(ibin).sumW() - background[itype][ipt]); - } - - // Integrate near-side yield - unsigned int lowerBin = _histYield[itype][ipt]->binIndexAt(-0.7 + 0.02); - unsigned int upperBin = _histYield[itype][ipt]->binIndexAt( 0.7 - 0.02) + 1; - nbins = upperBin - lowerBin; - numberOfBins[itype][ipt][0] = nbins; - nearSide[itype][ipt] = _histYield[itype][ipt]->integralRange(lowerBin, - upperBin) - nbins * background[itype][ipt]; - numberOfEntries[itype][ipt][0] = _histYield[itype][ipt]->integralRange( - lowerBin, upperBin) * _histTriggerCounter->bin(itype).sumW(); - - // Integrate away-side yield - lowerBin = _histYield[itype][ipt]->binIndexAt(M_PI - 0.7 + 0.02); - upperBin = _histYield[itype][ipt]->binIndexAt(M_PI + 0.7 - 0.02) + 1; - nbins = upperBin - lowerBin; - numberOfBins[itype][ipt][1] = nbins; - awaySide[itype][ipt] = _histYield[itype][ipt]->integralRange( - lowerBin, upperBin) - nbins * background[itype][ipt]; - numberOfEntries[itype][ipt][1] = - _histYield[itype][ipt]->integralRange(lowerBin, upperBin) * - _histTriggerCounter->bin(itype).sumW(); - - } - } - - // Variables for IAA/ICP plots - double dI = 0.0; - int near = 0; - int away = 1; - double xval[PT_BINS] = { 3.5, 5.0, 7.0, 9.0 }; - double xerr[PT_BINS] = { 0.5, 1.0, 1.0, 1.0 }; - - int types1[3] = {1, 2, 1}; - int types2[3] = {0, 0, 2}; - // Fill IAA/ICP plots for near side peak - for (int ihist = 0; ihist < 3; ihist++) { - int type1 = types1[ihist]; - int type2 = types2[ihist]; - for (int ipt = 0; ipt < PT_BINS; ipt++) { - dI = scalingFactor[type1] * scalingFactor[type1] * numberOfEntries[type1][ipt][near] + - scalingFactor[type2] * scalingFactor[type2] * numberOfEntries[type2][ipt][near] * - nearSide[type1][ipt] * nearSide[type1][ipt] / (nearSide[type2][ipt] * nearSide[type2][ipt]) + - numberOfBins[type1][ipt][near] * numberOfBins[type1][ipt][near] * backgroundError[type1][ipt] * - backgroundError[type1][ipt] + numberOfBins[type2][ipt][near] * numberOfBins[type2][ipt][near] * - backgroundError[type2][ipt] * backgroundError[type2][ipt] * nearSide[type1][ipt] * - nearSide[type1][ipt] / (nearSide[type2][ipt] * nearSide[type2][ipt]); + // Variables for near and away side peak calculation + double nearSide[EVENT_TYPES][PT_BINS] = { {0.0} }; + double awaySide[EVENT_TYPES][PT_BINS] = { {0.0} }; - dI = sqrt(dI)/nearSide[type2][ipt]; - _histIAA[ihist]->addPoint(xval[ipt], nearSide[type1][ipt] / nearSide[type2][ipt], xerr[ipt], dI); - } - } - - // Fill IAA/ICP plots for away side peak - for (int ihist = 0; ihist < 3; ihist++) { - int type1 = types1[ihist]; - int type2 = types2[ihist]; - for (int ipt = 0; ipt < PT_BINS; ipt++) { - dI = scalingFactor[type1] * scalingFactor[type1] * numberOfEntries[type1][ipt][away] + - scalingFactor[type2] * scalingFactor[type2] * numberOfEntries[type2][ipt][away] * - awaySide[type1][ipt] * awaySide[type1][ipt] / (awaySide[type2][ipt] * awaySide[type2][ipt]) + - numberOfBins[type1][ipt][away] * numberOfBins[type1][ipt][away] * backgroundError[type1][ipt] * - backgroundError[type1][ipt] + numberOfBins[type2][ipt][away] * numberOfBins[type2][ipt][away] * - backgroundError[type2][ipt] * backgroundError[type2][ipt] * awaySide[type1][ipt] * - awaySide[type1][ipt] / (awaySide[type2][ipt] * awaySide[type2][ipt]); + // Variables for background error calculation + double background[EVENT_TYPES][PT_BINS] = { {0.0} }; + double backgroundError[EVENT_TYPES][PT_BINS] = { {0.0} }; - dI = sqrt(dI)/awaySide[type2][ipt]; - _histIAA[ihist + 3]->addPoint(xval[ipt], awaySide[type1][ipt] / awaySide[type2][ipt], xerr[ipt], dI); - } - } + // Variables for integration error calculation + double scalingFactor[EVENT_TYPES] = {0.0}; + double numberOfEntries[EVENT_TYPES][PT_BINS][2] = { { {0.0} } }; + int numberOfBins[EVENT_TYPES][PT_BINS][2] = { { {0} } }; + + // For each event type + for (int itype = 0; itype < EVENT_TYPES; itype++) { + // For each pT range + for (int ipt = 0; ipt < PT_BINS; ipt++) { + // Check if histograms are fine + if (_histTriggerCounter->numEntries() == 0 || _histYield[itype][ipt]->numEntries() == 0) { + MSG_WARNING("There are no entries in one of the histograms"); + continue; + } + + // Scale yield histogram + if ((_histTriggerCounter->bin(itype).sumW() != 0)) { + scalingFactor[itype] = 1. / _histTriggerCounter->bin(itype).sumW(); + scale(_histYield[itype][ipt], (1. / _histTriggerCounter->bin(itype).sumW())); + } + + // Calculate background + double sum = 0.0; + int nbins = 0; + for (size_t ibin = 0; ibin < _histYield[itype][ipt]->numBins(); ibin++) { + if ((_histYield[itype][ipt]->bin(ibin).xMid() > (-0.5 * M_PI) && + _histYield[itype][ipt]->bin(ibin).xMid() < (-0.5 * M_PI + 0.4)) || + (_histYield[itype][ipt]->bin(ibin).xMid() > (0.5 * M_PI - 0.4) && + _histYield[itype][ipt]->bin(ibin).xMid() < (0.5 * M_PI + 0.4)) || + (_histYield[itype][ipt]->bin(ibin).xMid() > (1.5 * M_PI - 0.4) && + _histYield[itype][ipt]->bin(ibin).xMid() < (1.5 * M_PI))) { + sum += _histYield[itype][ipt]->bin(ibin).sumW(); + nbins += 1; + } + } + if (nbins == 0) { + MSG_WARNING("Failed to estimate background!"); + continue; + } + background[itype][ipt] = sum / nbins; + + // Calculate background error + sum = 0.0; + nbins = 0; + for (size_t ibin = 0; ibin < _histYield[itype][ipt]->numBins(); ibin++) { + if (_histYield[itype][ipt]->bin(ibin).xMid() > (0.5 * M_PI - 0.4) && + _histYield[itype][ipt]->bin(ibin).xMid() < (0.5 * M_PI + 0.4)) { + sum += (_histYield[itype][ipt]->bin(ibin).sumW() - background[itype][ipt]) * + (_histYield[itype][ipt]->bin(ibin).sumW() - background[itype][ipt]); + nbins++; + } + } + backgroundError[itype][ipt] = sqrt(sum / (nbins - 1)); + + // Fill histograms with removed background + for (size_t ibin = 0; ibin < _histYield[itype][ipt]->numBins(); ibin++) { + _histYieldBkgRemoved[itype][ipt]->fillBin(ibin, + _histYield[itype][ipt]->bin(ibin).sumW() - background[itype][ipt]); + } + + // Integrate near-side yield + size_t lowerBin = _histYield[itype][ipt]->binIndexAt(-0.7 + 0.02); + size_t upperBin = _histYield[itype][ipt]->binIndexAt( 0.7 - 0.02) + 1; + nbins = upperBin - lowerBin; + numberOfBins[itype][ipt][0] = nbins; + nearSide[itype][ipt] = _histYield[itype][ipt]->integralRange(lowerBin, upperBin) - nbins * background[itype][ipt]; + numberOfEntries[itype][ipt][0] = _histYield[itype][ipt]->integralRange(lowerBin, upperBin) * _histTriggerCounter->bin(itype).sumW(); + + // Integrate away-side yield + lowerBin = _histYield[itype][ipt]->binIndexAt(M_PI - 0.7 + 0.02); + upperBin = _histYield[itype][ipt]->binIndexAt(M_PI + 0.7 - 0.02) + 1; + nbins = upperBin - lowerBin; + numberOfBins[itype][ipt][1] = nbins; + awaySide[itype][ipt] = _histYield[itype][ipt]->integralRange(lowerBin, upperBin) - nbins * background[itype][ipt]; + numberOfEntries[itype][ipt][1] = _histYield[itype][ipt]->integralRange(lowerBin, upperBin) * _histTriggerCounter->bin(itype).sumW(); + + } + } + + // Variables for IAA/ICP plots + double dI = 0.0; + int near = 0; + int away = 1; + double xval[PT_BINS] = { 3.5, 5.0, 7.0, 9.0 }; + double xerr[PT_BINS] = { 0.5, 1.0, 1.0, 1.0 }; + + int types1[3] = {1, 2, 1}; + int types2[3] = {0, 0, 2}; + + // Fill IAA/ICP plots for near side peak + for (int ihist = 0; ihist < 3; ihist++) { + int type1 = types1[ihist]; + int type2 = types2[ihist]; + for (int ipt = 0; ipt < PT_BINS; ipt++) { + dI = scalingFactor[type1] * scalingFactor[type1] * numberOfEntries[type1][ipt][near] + + scalingFactor[type2] * scalingFactor[type2] * numberOfEntries[type2][ipt][near] * + nearSide[type1][ipt] * nearSide[type1][ipt] / (nearSide[type2][ipt] * nearSide[type2][ipt]) + + numberOfBins[type1][ipt][near] * numberOfBins[type1][ipt][near] * backgroundError[type1][ipt] * + backgroundError[type1][ipt] + numberOfBins[type2][ipt][near] * numberOfBins[type2][ipt][near] * + backgroundError[type2][ipt] * backgroundError[type2][ipt] * nearSide[type1][ipt] * + nearSide[type1][ipt] / (nearSide[type2][ipt] * nearSide[type2][ipt]); + + dI = sqrt(dI)/nearSide[type2][ipt]; + _histIAA[ihist]->addPoint(xval[ipt], nearSide[type1][ipt] / nearSide[type2][ipt], xerr[ipt], dI); + } + } + + // Fill IAA/ICP plots for away side peak + for (int ihist = 0; ihist < 3; ihist++) { + int type1 = types1[ihist]; + int type2 = types2[ihist]; + for (int ipt = 0; ipt < PT_BINS; ipt++) { + dI = scalingFactor[type1] * scalingFactor[type1] * numberOfEntries[type1][ipt][away] + + scalingFactor[type2] * scalingFactor[type2] * numberOfEntries[type2][ipt][away] * + awaySide[type1][ipt] * awaySide[type1][ipt] / (awaySide[type2][ipt] * awaySide[type2][ipt]) + + numberOfBins[type1][ipt][away] * numberOfBins[type1][ipt][away] * backgroundError[type1][ipt] * + backgroundError[type1][ipt] + numberOfBins[type2][ipt][away] * numberOfBins[type2][ipt][away] * + backgroundError[type2][ipt] * backgroundError[type2][ipt] * awaySide[type1][ipt] * + awaySide[type1][ipt] / (awaySide[type2][ipt] * awaySide[type2][ipt]); + + dI = sqrt(dI)/awaySide[type2][ipt]; + _histIAA[ihist + 3]->addPoint(xval[ipt], awaySide[type1][ipt] / awaySide[type2][ipt], xerr[ipt], dI); + } + } } - + } - + + private: - + static const int PT_BINS = 4; static const int EVENT_TYPES = 3; Histo1DPtr _histYield[EVENT_TYPES][PT_BINS]; Histo1DPtr _histYieldBkgRemoved[EVENT_TYPES][PT_BINS]; Histo1DPtr _histTriggerCounter; Scatter2DPtr _histIAA[6]; double pt_limits[5]; int event_type; string event_string[EVENT_TYPES + 1]; bool reentrant_flag; + }; - // The hook for the plugin system + DECLARE_RIVET_PLUGIN(ALICE_2012_I930312); + } diff --git a/analyses/pluginALICE/ALICE_2016_I1507157.cc b/analyses/pluginALICE/ALICE_2016_I1507157.cc --- a/analyses/pluginALICE/ALICE_2016_I1507157.cc +++ b/analyses/pluginALICE/ALICE_2016_I1507157.cc @@ -1,202 +1,167 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/AliceCommon.hh" #include "Rivet/Projections/PrimaryParticles.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/EventMixingFinalState.hh" namespace Rivet { - /// @brief Correlations of identified particles in pp. + /// @brief ALICE correlations of identified particles in pp /// Also showcasing use of EventMixingFinalState. - class ALICE_2016_I1507157 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ALICE_2016_I1507157); /// @name Analysis methods //@{ - - /// @brief Calculate angular distance between particles. - double phaseDif(double a1, double a2){ - double dif = a1 - a2; - while (dif < -M_PI/2) - dif += 2*M_PI; - while (dif > 3*M_PI/2) - dif -= 2*M_PI; - return dif; - } /// Book histograms and initialise projections before the run void init() { double etamax = 0.8; - double pTmin = 0.5; // GeV - + double pTmin = 0.5; // GeV + // Trigger declare(ALICE::V0AndTrigger(), "V0-AND"); // Charged tracks used to manage the mixing observable. ChargedFinalState cfsMult(Cuts::abseta < etamax); addProjection(cfsMult, "CFSMult"); - + // Primary particles. - PrimaryParticles pp({Rivet::PID::PIPLUS, Rivet::PID::KPLUS, - Rivet::PID::K0S, Rivet::PID::K0L, Rivet::PID::PROTON, - Rivet::PID::NEUTRON, Rivet::PID::LAMBDA, Rivet::PID::SIGMAMINUS, - Rivet::PID::SIGMAPLUS, Rivet::PID::XIMINUS, Rivet::PID::XI0, - Rivet::PID::OMEGAMINUS},Cuts::abseta < etamax && Cuts::pT > pTmin*GeV); + PrimaryParticles pp({Rivet::PID::PIPLUS, Rivet::PID::KPLUS, + Rivet::PID::K0S, Rivet::PID::K0L, Rivet::PID::PROTON, + Rivet::PID::NEUTRON, Rivet::PID::LAMBDA, Rivet::PID::SIGMAMINUS, + Rivet::PID::SIGMAPLUS, Rivet::PID::XIMINUS, Rivet::PID::XI0, + Rivet::PID::OMEGAMINUS},Cuts::abseta < etamax && Cuts::pT > pTmin*GeV); addProjection(pp,"APRIM"); // The event mixing projection declare(EventMixingFinalState(cfsMult, pp, 5, 0, 100, 10),"EVM"); // The particle pairs. pid = {{211, -211}, {321, -321}, {2212, -2212}, {3122, -3122}, {211, 211}, {321, 321}, {2212, 2212}, {3122, 3122}, {2212, 3122}, {2212, -3122}}; // The associated histograms in the data file. vector refdata = {"d04-x01-y01","d04-x01-y02","d04-x01-y03", - "d06-x01-y02","d05-x01-y01","d05-x01-y02","d05-x01-y03","d06-x01-y01", - "d01-x01-y02","d02-x01-y02"}; + "d06-x01-y02","d05-x01-y01","d05-x01-y02","d05-x01-y03","d06-x01-y01", + "d01-x01-y02","d02-x01-y02"}; for (int i = 0, N = refdata.size(); i < N; ++i) { // The ratio plots. - ratio.push_back(bookScatter2D(refdata[i], true)); - // Signal and mixed background. + ratio.push_back(bookScatter2D(refdata[i], true)); + // Signal and mixed background. signal.push_back(bookHisto1D("/TMP/" + refdata[i] + - "-s", *ratio[i], refdata[i] + "-s")); - background.push_back(bookHisto1D("/TMP/" + refdata[i] + - "-b", *ratio[i], refdata[i] + "-b")); + "-s", *ratio[i], refdata[i] + "-s")); + background.push_back(bookHisto1D("/TMP/" + refdata[i] + + "-b", *ratio[i], refdata[i] + "-b")); // Number of signal and mixed pairs. - nsp.push_back(0.); + nsp.push_back(0.); nmp.push_back(0.); - } + } } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = event.weight(); - + // Triggering if (!apply(event, "V0-AND")()) return; - // The projections - const PrimaryParticles& pp = applyProjection(event,"APRIM"); - const EventMixingFinalState& evm = - applyProjection(event, "EVM"); - // Get all mixing events - vector mixEvents = evm.getMixingEvents(); - - // If there are not enough events to mix, don't fill any histograms - if(mixEvents.size() == 0) - return; + // The primary particles + const PrimaryParticles& pp = apply(event, "APRIM"); + const Particles pparticles = pp.particles(); + + // The mixed events + const EventMixingFinalState& evm = apply(event, "EVM"); + const vector mixEvents = evm.getMixingEvents(); + if (mixEvents.empty()) vetoEvent; // Make a vector of mixed event particles vector mixParticles; size_t pSize = 0; - for(size_t i = 0; i < mixEvents.size(); ++i) - pSize+=mixEvents[i].size(); + for (size_t i = 0; i < mixEvents.size(); ++i) + pSize += mixEvents[i].size(); mixParticles.reserve(pSize); - for(size_t i = 0; i < mixEvents.size(); ++i) - mixParticles.insert(mixParticles.end(), mixEvents[i].begin(), - mixEvents[i].end()); - - // Shuffle the particles in the mixing events + for (size_t i = 0; i < mixEvents.size(); ++i) + mixParticles.insert(mixParticles.end(), mixEvents[i].begin(), mixEvents[i].end()); random_shuffle(mixParticles.begin(), mixParticles.end()); - - for(const Particle& p1 : pp.particles()) { + + for (size_t ip1 = 0; ip1 < pparticles.size()-1; ++ip1) { + const Particle& p1 = pparticles[ip1]; + // Start by doing the signal distributions - for(const Particle& p2 : pp.particles()) { - if(p1 == p2) - continue; - double dEta = abs(p1.eta() - p2.eta()); - double dPhi = phaseDif(p1.phi(), p2.phi()); - if(dEta < 1.3) { - for (int i = 0, N = pid.size(); i < N; ++i) { - int pid1 = pid[i].first; - int pid2 = pid[i].second; - bool samesign = (pid1 * pid2 > 0); - if (samesign && ((pid1 == p1.pid() && pid2 == p2.pid()) || - (pid1 == -p1.pid() && pid2 == -p2.pid()))) { - signal[i]->fill(dPhi, weight); - nsp[i] += 1.0; - } - if (!samesign && abs(pid1) == abs(pid2) && - pid1 == p1.pid() && pid2 == p2.pid()) { - signal[i]->fill(dPhi, weight); - nsp[i] += 1.0; - } - if (!samesign && abs(pid1) != abs(pid2) && - ( (pid1 == p1.pid() && pid2 == p2.pid()) || - (pid2 == p1.pid() && pid1 == p2.pid()) ) ) { - signal[i]->fill(dPhi, weight); - nsp[i] += 1.0; - } - } - } - } - // Then do the background distribution - for(const Particle& pMix : mixParticles){ - double dEta = abs(p1.eta() - pMix.eta()); - double dPhi = phaseDif(p1.phi(), pMix.phi()); - if(dEta < 1.3) { - for (int i = 0, N = pid.size(); i < N; ++i) { - int pid1 = pid[i].first; - int pid2 = pid[i].second; - bool samesign = (pid1 * pid2 > 0); - if (samesign && ((pid1 == p1.pid() && pid2 == pMix.pid()) || - (pid1 == -p1.pid() && pid2 == -pMix.pid()))) { - background[i]->fill(dPhi, weight); - nmp[i] += 1.0; - } - if (!samesign && abs(pid1) == abs(pid2) && - pid1 == p1.pid() && pid2 == pMix.pid()) { - background[i]->fill(dPhi, weight); - nmp[i] += 1.0; - } - if (!samesign && abs(pid1) != abs(pid2) && - ( (pid1 == p1.pid() && pid2 == pMix.pid()) || - (pid2 == p1.pid() && pid1 == pMix.pid()) ) ) { - background[i]->fill(dPhi, weight); - nmp[i] += 1.0; - } - } - } - } + for (size_t ip2 = 0; ip2 < pparticles.size(); ++ip1) { + if (ip1 == ip2) continue; + const Particle& p2 = pparticles[ip2]; + const double dEta = deltaEta(p1, p2); + const double dPhi = deltaPhi(p1, p2, true); + if (dEta > 1.3) continue; + for (int i = 0, N = pid.size(); i < N; ++i) { + const int pid1 = pid[i].first; + const int pid2 = pid[i].second; + const bool samesign = (pid1 * pid2 > 0); + const bool pidmatch1 = (pid1 == p1.pid() && pid2 == p2.pid()) || (pid1 == -p1.pid() && pid2 == -p2.pid()); + const bool pidmatch2 = abs(pid1) == abs(pid2) && pid1 == p1.pid() && pid2 == p2.pid(); + const bool pidmatch3 = abs(pid1) != abs(pid2) && ( (pid1 == p1.pid() && pid2 == p2.pid()) || (pid2 == p1.pid() && pid1 == p2.pid()) ); + if ((samesign && pidmatch1) || (!samesign && (pidmatch2 || pidmatch3))) { + signal[i]->fill(dPhi, weight); + nsp[i] += 1.0; + } + } + } + + // Then do the background distribution + for (const Particle& pMix : mixParticles){ + const double dEta = deltaEta(p1, pMix); + const double dPhi = deltaPhi(p1, pMix, true); + if (dEta > 1.3) continue; + for (int i = 0, N = pid.size(); i < N; ++i) { + const int pid1 = pid[i].first; + const int pid2 = pid[i].second; + const bool samesign = (pid1 * pid2 > 0); + const bool pidmatch1 = (pid1 == p1.pid() && pid2 == pMix.pid()) || (pid1 == -p1.pid() && pid2 == -pMix.pid()); + const bool pidmatch2 = abs(pid1) == abs(pid2) && pid1 == p1.pid() && pid2 == pMix.pid(); + const bool pidmatch3 = abs(pid1) != abs(pid2) && ( (pid1 == p1.pid() && pid2 == pMix.pid()) || (pid2 == p1.pid() && pid1 == pMix.pid()) ); + if ((samesign && pidmatch1) || (!samesign && (pidmatch2 || pidmatch3))) { + background[i]->fill(dPhi, weight); + nmp[i] += 1.0; + } + } + } } } /// Normalise histograms etc., after the run void finalize() { for (int i = 0, N = pid.size(); i < N; ++i) { - double sc = nmp[i] / nsp[i]; - signal[i]->scaleW(sc); - divide(signal[i],background[i],ratio[i]); + const double sc = nmp[i] / nsp[i]; + signal[i]->scaleW(sc); + divide(signal[i],background[i],ratio[i]); } } //@} /// @name Histograms //@{ vector > pid; vector signal; vector background; vector ratio; vector nsp; vector nmp; //@} }; - // The hook for the plugin system DECLARE_RIVET_PLUGIN(ALICE_2016_I1507157); - } diff --git a/analyses/pluginATLAS/ATLAS_2012_I1203852.cc b/analyses/pluginATLAS/ATLAS_2012_I1203852.cc --- a/analyses/pluginATLAS/ATLAS_2012_I1203852.cc +++ b/analyses/pluginATLAS/ATLAS_2012_I1203852.cc @@ -1,373 +1,375 @@ // -*- 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 { + /// 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(); } }; - /// @name ZZ analysis + /// ZZ analysis class ATLAS_2012_I1203852 : public Analysis { public: /// Default constructor ATLAS_2012_I1203852() : Analysis("ATLAS_2012_I1203852") { } + void init() { // NB Missing ET is not required to be neutrinos FinalState fs(-5.0, 5.0, 0.0*GeV); // Final states to form Z bosons vids.push_back(make_pair(PID::ELECTRON, PID::POSITRON)); vids.push_back(make_pair(PID::MUON, PID::ANTIMUON)); IdentifiedFinalState Photon(fs); Photon.acceptIdPair(PID::PHOTON); IdentifiedFinalState bare_EL(fs); bare_EL.acceptIdPair(PID::ELECTRON); IdentifiedFinalState bare_MU(fs); bare_MU.acceptIdPair(PID::MUON); // Selection 1: ZZ-> llll selection Cut etaranges_lep = Cuts::abseta < 3.16 && Cuts::pT > 7*GeV; DressedLeptons electron_sel4l(Photon, bare_EL, 0.1, etaranges_lep); declare(electron_sel4l, "ELECTRON_sel4l"); DressedLeptons muon_sel4l(Photon, bare_MU, 0.1, etaranges_lep); declare(muon_sel4l, "MUON_sel4l"); // Selection 2: ZZ-> llnunu selection Cut etaranges_lep2 = Cuts::abseta < 2.5 && Cuts::pT > 10*GeV; DressedLeptons electron_sel2l2nu(Photon, bare_EL, 0.1, etaranges_lep2); declare(electron_sel2l2nu, "ELECTRON_sel2l2nu"); DressedLeptons muon_sel2l2nu(Photon, bare_MU, 0.1, etaranges_lep2); declare(muon_sel2l2nu, "MUON_sel2l2nu"); /// Get all neutrinos. These will not be used to form jets. IdentifiedFinalState neutrino_fs(Cuts::abseta < 4.5); neutrino_fs.acceptNeutrinos(); declare(neutrino_fs, "NEUTRINO_FS"); // Calculate missing ET from the visible final state, not by requiring neutrinos addProjection(MissingMomentum(Cuts::abseta < 4.5), "MISSING"); VetoedFinalState jetinput; jetinput.addVetoOnThisFinalState(bare_MU); jetinput.addVetoOnThisFinalState(neutrino_fs); FastJets jetpro(fs, FastJets::ANTIKT, 0.4); declare(jetpro, "jet"); // Both ZZ on-shell histos _h_ZZ_xsect = bookHisto1D(1, 1, 1); _h_ZZ_ZpT = bookHisto1D(3, 1, 1); _h_ZZ_phill = bookHisto1D(5, 1, 1); _h_ZZ_mZZ = bookHisto1D(7, 1, 1); // One Z off-shell (ZZstar) histos _h_ZZs_xsect = bookHisto1D(1, 1, 2); // ZZ -> llnunu histos _h_ZZnunu_xsect = bookHisto1D(1, 1, 3); _h_ZZnunu_ZpT = bookHisto1D(4, 1, 1); _h_ZZnunu_phill = bookHisto1D(6, 1, 1); _h_ZZnunu_mZZ = bookHisto1D(8, 1, 1); } /// Do the analysis void analyze(const Event& e) { const double weight = e.weight(); //////////////////////////////////////////////////////////////////// // 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(); 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() ); //////////////////////////////////////////////////////////////////// // OVERLAP removal dR(l,l)>0.2 //////////////////////////////////////////////////////////////////// - foreach ( const DressedLepton& l1, leptonsFS_sel4l) { + for (const DressedLepton& l1 : leptonsFS_sel4l) { bool isolated = true; - foreach (DressedLepton& l2, leptonsFS_sel4l) { + for (DressedLepton& l2 : leptonsFS_sel4l) { const double dR = deltaR(l1, l2); - if (dR < 0.2 && l1 != l2) { isolated = false; break; } + if (dR < 0.2 && !isSame(l1, l2)) { isolated = false; break; } } if (isolated) leptons_sel4l.push_back(l1); } ////////////////////////////////////////////////////////////////// // Exactly two opposite charged leptons ////////////////////////////////////////////////////////////////// // calculate total 'flavour' charge double totalcharge = 0; - foreach (Particle& l, leptons_sel4l) totalcharge += l.pid(); + for (const Particle& l : leptons_sel4l) totalcharge += l.pid(); // Analyze 4 lepton events if (leptons_sel4l.size() == 4 && totalcharge == 0 ) { Zstate Z1, Z2; // Identifies Z states from 4 lepton pairs identifyZstates(Z1, Z2,leptons_sel4l); //////////////////////////////////////////////////////////////////////////// // Z MASS WINDOW // -ZZ: for both Z: 6620 GeV /////////////////////////////////////////////////////////////////////////// Zstate leadPtZ = std::max(Z1, Z2, Zstate::cmppT); double mZ1 = Z1.mom().mass(); double mZ2 = Z2.mom().mass(); double ZpT = leadPtZ.mom().pT(); double phill = fabs(deltaPhi(leadPtZ.first, leadPtZ.second)); if (phill > M_PI) phill = 2*M_PI-phill; double mZZ = (Z1.mom() + Z2.mom()).mass(); if (mZ1 > 20*GeV && mZ2 > 20*GeV) { // ZZ* selection if (inRange(mZ1, 66*GeV, 116*GeV) || inRange(mZ2, 66*GeV, 116*GeV)) { _h_ZZs_xsect -> fill(sqrtS()*GeV, weight); } // ZZ selection if (inRange(mZ1, 66*GeV, 116*GeV) && inRange(mZ2, 66*GeV, 116*GeV)) { _h_ZZ_xsect -> fill(sqrtS()*GeV, weight); _h_ZZ_ZpT -> fill(ZpT , weight); _h_ZZ_phill -> fill(phill , weight); _h_ZZ_mZZ -> fill(mZZ , weight); } } } //////////////////////////////////////////////////////////////////// /// preselection of leptons for ZZ-> llnunu final state //////////////////////////////////////////////////////////////////// Particles leptons_sel2l2nu; // output const vector& mu_sel2l2nu = apply(e, "MUON_sel2l2nu").dressedLeptons(); const vector& el_sel2l2nu = apply(e, "ELECTRON_sel2l2nu").dressedLeptons(); vector leptonsFS_sel2l2nu; leptonsFS_sel2l2nu.insert( leptonsFS_sel2l2nu.end(), mu_sel2l2nu.begin(), mu_sel2l2nu.end() ); leptonsFS_sel2l2nu.insert( leptonsFS_sel2l2nu.end(), el_sel2l2nu.begin(), el_sel2l2nu.end() ); // Lepton preselection for ZZ-> llnunu if ((mu_sel2l2nu.empty() || el_sel2l2nu.empty()) // cannot have opposite flavour && (leptonsFS_sel2l2nu.size() == 2) // exactly two leptons && (leptonsFS_sel2l2nu[0].charge() * leptonsFS_sel2l2nu[1].charge() < 1 ) // opposite charge && (deltaR(leptonsFS_sel2l2nu[0], leptonsFS_sel2l2nu[1]) > 0.3) // overlap removal && (leptonsFS_sel2l2nu[0].pT() > 20*GeV && leptonsFS_sel2l2nu[1].pT() > 20*GeV)) { // trigger requirement leptons_sel2l2nu.insert(leptons_sel2l2nu.end(), leptonsFS_sel2l2nu.begin(), leptonsFS_sel2l2nu.end()); } if (leptons_sel2l2nu.empty()) vetoEvent; // no further analysis, fine to veto Particles leptons_sel2l2nu_jetveto; foreach (const DressedLepton& l, mu_sel2l2nu) leptons_sel2l2nu_jetveto.push_back(l.constituentLepton()); foreach (const DressedLepton& l, el_sel2l2nu) leptons_sel2l2nu_jetveto.push_back(l.constituentLepton()); double ptll = (leptons_sel2l2nu[0].momentum() + leptons_sel2l2nu[1].momentum()).pT(); // Find Z1-> ll FinalState fs2(-3.2, 3.2); InvMassFinalState imfs(fs2, vids, 20*GeV, sqrtS()); imfs.calc(leptons_sel2l2nu); if (imfs.particlePairs().size() != 1) vetoEvent; const ParticlePair& Z1constituents = imfs.particlePairs()[0]; FourMomentum Z1 = Z1constituents.first.momentum() + Z1constituents.second.momentum(); // Z to neutrinos candidate from missing ET const MissingMomentum & missmom = applyProjection(e, "MISSING"); const FourMomentum Z2 = missmom.missingMomentum(ZMASS); double met_Znunu = missmom.missingEt(); //Z2.pT(); // mTZZ const double mT2_1st_term = add_quad(ZMASS, ptll) + add_quad(ZMASS, met_Znunu); const double mT2_2nd_term = Z1.px() + Z2.px(); const double mT2_3rd_term = Z1.py() + Z2.py(); const double mTZZ = sqrt(sqr(mT2_1st_term) - sqr(mT2_2nd_term) - sqr(mT2_3rd_term)); if (!inRange(Z2.mass(), 66*GeV, 116*GeV)) vetoEvent; if (!inRange(Z1.mass(), 76*GeV, 106*GeV)) vetoEvent; ///////////////////////////////////////////////////////////// // AXIAL MET < 75 GeV //////////////////////////////////////////////////////////// double dPhiZ1Z2 = fabs(deltaPhi(Z1, Z2)); if (dPhiZ1Z2 > M_PI) dPhiZ1Z2 = 2*M_PI - dPhiZ1Z2; const double axialEtmiss = -Z2.pT()*cos(dPhiZ1Z2); if (axialEtmiss < 75*GeV) vetoEvent; const double ZpT = Z1.pT(); double phill = fabs(deltaPhi(Z1constituents.first, Z1constituents.second)); if (phill > M_PI) phill = 2*M_PI - phill; //////////////////////////////////////////////////////////////////////////// // JETS // -"j": found by "jetpro" projection && pT() > 25 GeV && |eta| < 4.5 // -"goodjets": "j" && dR(electron/muon,jet) > 0.3 // // JETVETO: veto all events with at least one good jet /////////////////////////////////////////////////////////////////////////// vector good_jets; foreach (const Jet& j, apply(e, "jet").jetsByPt(25)) { if (j.abseta() > 4.5) continue; bool isLepton = 0; foreach (const Particle& l, leptons_sel2l2nu_jetveto) { const double dR = deltaR(l.momentum(), j.momentum()); if (dR < 0.3) { isLepton = true; break; } } if (!isLepton) good_jets.push_back(j); } size_t n_sel_jets = good_jets.size(); if (n_sel_jets != 0) vetoEvent; ///////////////////////////////////////////////////////////// // Fractional MET and lepton pair difference: "RatioMet"< 0.4 //////////////////////////////////////////////////////////// double ratioMet = fabs(Z2.pT() - Z1.pT()) / Z1.pT(); if (ratioMet > 0.4 ) vetoEvent; // End of ZZllnunu selection: now fill histograms _h_ZZnunu_xsect->fill(sqrtS()/GeV, weight); _h_ZZnunu_ZpT ->fill(ZpT, weight); _h_ZZnunu_phill->fill(phill, weight); _h_ZZnunu_mZZ ->fill(mTZZ, weight); } /// Finalize void finalize() { const double norm = crossSection()/sumOfWeights()/femtobarn; scale(_h_ZZ_xsect, norm); normalize(_h_ZZ_ZpT); normalize(_h_ZZ_phill); normalize(_h_ZZ_mZZ); scale(_h_ZZs_xsect, norm); scale(_h_ZZnunu_xsect, norm); normalize(_h_ZZnunu_ZpT); normalize(_h_ZZnunu_phill); normalize(_h_ZZnunu_mZZ); } private: void identifyZstates(Zstate& Z1, Zstate& Z2, const Particles& leptons_sel4l); Histo1DPtr _h_ZZ_xsect, _h_ZZ_ZpT, _h_ZZ_phill, _h_ZZ_mZZ; Histo1DPtr _h_ZZs_xsect; Histo1DPtr _h_ZZnunu_xsect, _h_ZZnunu_ZpT, _h_ZZnunu_phill, _h_ZZnunu_mZZ; vector< pair > vids; const double ZMASS = 91.1876; // GeV }; /// 4l to ZZ assignment -- algorithm void ATLAS_2012_I1203852::identifyZstates(Zstate& Z1, Zstate& Z2, const Particles& leptons_sel4l) { ///////////////////////////////////////////////////////////////////////////// /// ZZ->4l pairing /// - Exactly two same flavour opposite charged leptons /// - Ambiguities in pairing are resolved by choosing the combination /// that results in the smaller value of the sum |mll - mZ| for the two pairs ///////////////////////////////////////////////////////////////////////////// Particles part_pos_el, part_neg_el, part_pos_mu, part_neg_mu; foreach (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); } } // ee/mm channel if ( part_neg_el.size() == 2 || part_neg_mu.size() == 2) { Zstate Zcand_1, Zcand_2, Zcand_3, Zcand_4; if (part_neg_el.size() == 2) { // ee Zcand_1 = Zstate( ParticlePair( part_neg_el[0], part_pos_el[0] ) ); Zcand_2 = Zstate( ParticlePair( part_neg_el[0], part_pos_el[1] ) ); Zcand_3 = Zstate( ParticlePair( part_neg_el[1], part_pos_el[0] ) ); Zcand_4 = Zstate( ParticlePair( part_neg_el[1], part_pos_el[1] ) ); } else { // mumu Zcand_1 = Zstate( ParticlePair( part_neg_mu[0], part_pos_mu[0] ) ); Zcand_2 = Zstate( ParticlePair( part_neg_mu[0], part_pos_mu[1] ) ); Zcand_3 = Zstate( ParticlePair( part_neg_mu[1], part_pos_mu[0] ) ); Zcand_4 = Zstate( ParticlePair( part_neg_mu[1], part_pos_mu[1] ) ); } // We can have the following pairs: (Z1 + Z4) || (Z2 + Z3) double minValue_1, minValue_2; minValue_1 = fabs( Zcand_1.mom().mass() - ZMASS ) + fabs( Zcand_4.mom().mass() - ZMASS); minValue_2 = fabs( Zcand_2.mom().mass() - ZMASS ) + fabs( Zcand_3.mom().mass() - ZMASS); if (minValue_1 < minValue_2 ) { Z1 = Zcand_1; Z2 = Zcand_4; } else { Z1 = Zcand_2; Z2 = Zcand_3; } // emu channel } else if (part_neg_mu.size() == 1 && part_neg_el.size() == 1) { Z1 = Zstate ( ParticlePair (part_neg_mu[0], part_pos_mu[0] ) ); Z2 = Zstate ( ParticlePair (part_neg_el[0], part_pos_el[0] ) ); } } // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2012_I1203852); } diff --git a/analyses/pluginATLAS/ATLAS_2013_I1190187.cc b/analyses/pluginATLAS/ATLAS_2013_I1190187.cc --- a/analyses/pluginATLAS/ATLAS_2013_I1190187.cc +++ b/analyses/pluginATLAS/ATLAS_2013_I1190187.cc @@ -1,261 +1,261 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/MissingMomentum.hh" namespace Rivet { /// ATLAS Wee Wemu Wmumu analysis at Z TeV class ATLAS_2013_I1190187 : public Analysis { public: /// Default constructor ATLAS_2013_I1190187() : Analysis("ATLAS_2013_I1190187") { } void init() { FinalState fs; Cut etaRanges_EL = (Cuts::abseta < 1.37 || Cuts::absetaIn(1.52, 2.47)) && Cuts::pT > 20*GeV; Cut etaRanges_MU = Cuts::abseta < 2.4 && Cuts::pT > 20*GeV; MissingMomentum met(fs); declare(met, "MET"); IdentifiedFinalState Photon(fs); Photon.acceptIdPair(PID::PHOTON); IdentifiedFinalState bare_EL(fs); bare_EL.acceptIdPair(PID::ELECTRON); IdentifiedFinalState bare_MU(fs); bare_MU.acceptIdPair(PID::MUON); IdentifiedFinalState neutrinoFS(fs); neutrinoFS.acceptNeutrinos(); declare(neutrinoFS, "Neutrinos"); //////////////////////////////////////////////////////// // DRESSED LEPTONS // 3.arg: 0.1 = dR(lep,phot) // 4.arg: true = do clustering // 7.arg: false = ignore photons from hadron or tau // ////////////////////////////////////////////////////////// DressedLeptons electronFS(Photon, bare_EL, 0.1, etaRanges_EL); declare(electronFS, "ELECTRON_FS"); DressedLeptons muonFS(Photon, bare_MU, 0.1, etaRanges_MU); declare(muonFS, "MUON_FS"); VetoedFinalState jetinput; jetinput.addVetoOnThisFinalState(bare_MU); jetinput.addVetoOnThisFinalState(neutrinoFS); FastJets jetpro(jetinput, FastJets::ANTIKT, 0.4); declare(jetpro, "jet"); // Book histograms _h_Wl1_pT_mumu = bookHisto1D(1, 1, 2); _h_Wl1_pT_ee = bookHisto1D(1, 1, 1); _h_Wl1_pT_emu = bookHisto1D(1, 1, 3); _h_Wl1_pT_inclusive = bookHisto1D(4, 1, 1); } /// Do the analysis void analyze(const Event& e) { const vector& muonFS = apply(e, "MUON_FS").dressedLeptons(); const vector& electronFS = apply(e, "ELECTRON_FS").dressedLeptons(); const MissingMomentum& met = apply(e, "MET"); vector dressed_lepton, isolated_lepton, fiducial_lepton; dressed_lepton.insert(dressed_lepton.end(), muonFS.begin(), muonFS.end()); dressed_lepton.insert(dressed_lepton.end(), electronFS.begin(), electronFS.end()); //////////////////////////////////////////////////////////////////////////// // OVERLAP REMOVAL // -electrons with dR(e,mu)<0.1 are removed // -lower pT electrons with dR(e,e)<0.1 are removed // //////////////////////////////////////////////////////////////////////////// - foreach (DressedLepton& l1, dressed_lepton) { + for (DressedLepton& l1 : dressed_lepton) { bool l_isolated = true; - foreach (DressedLepton& l2, dressed_lepton) { - if (l1 != l2 && l2.constituentLepton().abspid() == PID::ELECTRON) { + for (DressedLepton& l2 : dressed_lepton) { + if (!isSame(l1, l2) && l2.constituentLepton().abspid() == PID::ELECTRON) { double overlapControl_ll= deltaR(l1.constituentLepton(),l2.constituentLepton()); if (overlapControl_ll < 0.1) { l_isolated = false; // e/e overlap removal if (l1.constituentLepton().abspid() == PID::ELECTRON) { if (l1.constituentLepton().pT()>l2.constituentLepton().pT()) { isolated_lepton.push_back(l1);//keep e with highest pT } else { isolated_lepton.push_back(l2);//keep e with highest pT } } // e/mu overlap removal if (l1.constituentLepton().abspid() == PID::MUON) isolated_lepton.push_back(l1); //keep mu } } } if (l_isolated) isolated_lepton.push_back(l1); } /////////////////////////////////////////////////////////////////////////////////////////////////////////// // PRESELECTION: // "isolated_lepton:" // * electron: pt>20 GeV, |eta|<1.37, 1.52<|eta|<2.47, dR(electron,muon)>0.1 // * muon: pt>20 GeV, |eta|<2.4 // * dR(l,l)>0.1 // // "fiducial_lepton"= isolated_lepton with // * 2 leptons (e or mu) // * leading lepton pt (pT_l1) >25 GeV // * opposite charged leptons // /////////////////////////////////////////////////////////////////////////////////////////////////////////// if (isolated_lepton.size() != 2) vetoEvent; sort(isolated_lepton.begin(), isolated_lepton.end(), cmpMomByPt); if (isolated_lepton[0].pT() > 25*GeV && threeCharge(isolated_lepton[0]) != threeCharge(isolated_lepton[1])) { fiducial_lepton.insert(fiducial_lepton.end(), isolated_lepton.begin(), isolated_lepton.end()); } if (fiducial_lepton.size() == 0) vetoEvent; double pT_l1 = fiducial_lepton[0].pT(); double M_l1l2 = (fiducial_lepton[0].momentum() + fiducial_lepton[1].momentum()).mass(); double pT_l1l2 = (fiducial_lepton[0].momentum() + fiducial_lepton[1].momentum()).pT(); ///////////////////////////////////////////////////////////////////////// // JETS // -"alljets": found by "jetpro" projection && pT()>25 GeV && |y|<4.5 // -"vetojets": "alljets" && dR(electron,jet)>0.3 // ///////////////////////////////////////////////////////////////////////// Jets alljets, vetojets; foreach (const Jet& j, apply(e, "jet").jetsByPt(25)) { if (j.absrap() > 4.5 ) continue; alljets.push_back(j); bool deltaRcontrol = true; foreach (DressedLepton& fl,fiducial_lepton) { if (fl.constituentLepton().abspid() == PID::ELECTRON) { //electrons double deltaRjets = deltaR(fl.constituentLepton().momentum(), j.momentum(), RAPIDITY); if (deltaRjets <= 0.3) deltaRcontrol = false; //false if at least one electron is in the overlap region } } if (deltaRcontrol) vetojets.push_back(j); } ///////////////////////////////////////////////////////////////////////////////////////////////// // MISSING ETrel // -"mismom": fourvector of invisible momentum found by "met" projection // -"delta_phi": delta phi between mismom and the nearest "fiducial_lepton" or "vetojet" // -"MET_rel": missing transverse energy defined as: // *"mismom" for "delta_phi" >= (0.5*pi) // *"mismom.pT()*sin(delta_phi)" for "delta_phi" < (0.5*pi) // ///////////////////////////////////////////////////////////////////////////////////////////////// FourMomentum mismom; double MET_rel = 0, delta_phi = 0; mismom = -met.visibleMomentum(); vector vL_MET_angle, vJet_MET_angle; vL_MET_angle.push_back(fabs(deltaPhi(fiducial_lepton[0].momentum(), mismom))); vL_MET_angle.push_back(fabs(deltaPhi(fiducial_lepton[1].momentum(), mismom))); foreach (double& lM, vL_MET_angle) if (lM > M_PI) lM = 2*M_PI - lM; std::sort(vL_MET_angle.begin(), vL_MET_angle.end()); if (vetojets.size() == 0) delta_phi = vL_MET_angle[0]; if (vetojets.size() > 0) { foreach (Jet& vj, vetojets) { double jet_MET_angle = fabs(deltaPhi(vj.momentum(), mismom)); if (jet_MET_angle > M_PI) jet_MET_angle = 2*M_PI - jet_MET_angle; vJet_MET_angle.push_back(jet_MET_angle); } std::sort(vJet_MET_angle.begin(), vJet_MET_angle.end()); if (vL_MET_angle[0] <= vJet_MET_angle[0]) delta_phi = vL_MET_angle[0]; if (vL_MET_angle[0] > vJet_MET_angle[0]) delta_phi = vJet_MET_angle[0]; } if (delta_phi >= (0.5*M_PI)) delta_phi = 0.5*M_PI; MET_rel = mismom.pT()*sin(delta_phi); /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // CUTS // -jetveto: event with at least one vetojet is vetoed // -M_Z: Z mass M_Z=91.1876*GeV // // * ee channel: MET_rel > 45 GeV, M_l1l2 > 15 GeV, |M_l1l2-M_Z| > 15 GeV, jetveto, pT_l1l2 > 30 GeV // * mumu channel: MET_rel > 45 GeV, M_l1l2 > 15 GeV, |M_l1l2-M_Z| > 15 GeV, jetveto, pT_l1l2 > 30 GeV // * emu channel: MET_rel > 25 GeV, M_l1l2 > 10 GeV, |M_l1l2-M_Z| > 0 GeV, jetveto, pT_l1l2 > 30 GeV // //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // Get event weight for histo filling const double weight = e.weight(); // ee channel if (fiducial_lepton[0].abspid() == PID::ELECTRON && fiducial_lepton[1].abspid() == PID::ELECTRON) { if (MET_rel <= 45*GeV) vetoEvent; if (M_l1l2 <= 15*GeV) vetoEvent; if (fabs(M_l1l2 - 91.1876*GeV) <= 15*GeV) vetoEvent; if (vetojets.size() != 0) vetoEvent; if (pT_l1l2 <= 30*GeV) vetoEvent; _h_Wl1_pT_ee->fill(sqrtS()*GeV, weight); _h_Wl1_pT_inclusive->fill(pT_l1, weight); } // mumu channel else if (fiducial_lepton[0].abspid() == PID::MUON && fiducial_lepton[1].abspid() == PID::MUON) { if (MET_rel <= 45*GeV) vetoEvent; if (M_l1l2 <= 15*GeV) vetoEvent; if (fabs(M_l1l2-91.1876*GeV) <= 15*GeV) vetoEvent; if (vetojets.size() != 0) vetoEvent; if (pT_l1l2 <= 30*GeV) vetoEvent; _h_Wl1_pT_mumu->fill(sqrtS()*GeV, weight); _h_Wl1_pT_inclusive->fill(pT_l1, weight); } // emu channel else if (fiducial_lepton[0].abspid() != fiducial_lepton[1].abspid()) { if (MET_rel <= 25*GeV) vetoEvent; if (M_l1l2 <= 10*GeV) vetoEvent; if (vetojets.size() != 0) vetoEvent; if (pT_l1l2 <= 30*GeV) vetoEvent; _h_Wl1_pT_emu->fill(sqrtS()*GeV, weight); _h_Wl1_pT_inclusive->fill(pT_l1, weight); } } /// Finalize void finalize() { const double norm = crossSection()/sumOfWeights()/femtobarn; scale(_h_Wl1_pT_ee, norm); scale(_h_Wl1_pT_mumu, norm); scale(_h_Wl1_pT_emu, norm); normalize(_h_Wl1_pT_inclusive, 1); } private: Histo1DPtr _h_Wl1_pT_ee, _h_Wl1_pT_mumu, _h_Wl1_pT_emu, _h_Wl1_pT_inclusive; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2013_I1190187); } diff --git a/analyses/pluginATLAS/ATLAS_2017_I1625109.cc b/analyses/pluginATLAS/ATLAS_2017_I1625109.cc --- a/analyses/pluginATLAS/ATLAS_2017_I1625109.cc +++ b/analyses/pluginATLAS/ATLAS_2017_I1625109.cc @@ -1,309 +1,311 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { - + + class ATLAS_2017_I1625109 : public Analysis { public: /// Constructor /// @brief measurement of on-shell ZZ at 13 TeV DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2017_I1625109); /// @name Analysis methods //@{ struct Dilepton { Dilepton() {}; Dilepton(const ParticlePair & _leptons) : leptons(_leptons) {} FourMomentum momentum() const { return leptons.first.mom() + leptons.second.mom(); } ParticlePair leptons; }; - - + + struct Quadruplet { - + vector getLeptonsSortedByPt() const { - vector out = { leadingDilepton.leptons.first, leadingDilepton.leptons.second, + vector out = { leadingDilepton.leptons.first, leadingDilepton.leptons.second, subleadingDilepton.leptons.first, subleadingDilepton.leptons.second }; std::sort(out.begin(), out.end(), cmpMomByPt); return out; } - + Quadruplet(const Dilepton& dilepton1, const Dilepton& dilepton2) { if (dilepton1.momentum().pt() > dilepton2.momentum().pt()) { leadingDilepton = dilepton1; subleadingDilepton = dilepton2; } else { leadingDilepton = dilepton2; subleadingDilepton = dilepton1; } leptonsSortedByPt = getLeptonsSortedByPt(); } FourMomentum momentum() const { return leadingDilepton.momentum() + subleadingDilepton.momentum(); } - + double distanceFromZMass() const { return abs(leadingDilepton.momentum().mass() - Z_mass) + abs(subleadingDilepton.momentum().mass() - Z_mass); - } + } Dilepton leadingDilepton; Dilepton subleadingDilepton; vector leptonsSortedByPt; }; - + typedef vector Quadruplets; - + typedef std::pair IndexPair; - + + vector getOppositeChargePairsIndices(const vector& leptons) { vector indices = {}; if (leptons.size() < 2) return indices; for (size_t i = 0; i < leptons.size(); ++i) { for (size_t k = i+1; k < leptons.size(); ++k) { const auto charge_i = leptons.at(i).charge(); const auto charge_k = leptons.at(k).charge(); if (charge_i == -charge_k) { indices.push_back(std::make_pair(i, k)); } } } return indices; } - + bool indicesOverlap(IndexPair a, IndexPair b) { return (a.first == b.first || a.first == b.second || a.second == b.first || a.second == b.second); } - - + + bool passesHierarchicalPtRequirements(const Quadruplet& quadruplet) { const auto& sorted_leptons = quadruplet.leptonsSortedByPt; if (sorted_leptons.at(0).pt() < 20*GeV) return false; if (sorted_leptons.at(1).pt() < 15*GeV) return false; if (sorted_leptons.at(2).pt() < 10*GeV) return false; return true; } - + bool passesDileptonMinimumMassRequirement(const Quadruplet& quadruplet) { const auto& leptons = quadruplet.leptonsSortedByPt; - for (const auto& l1 : leptons) { - for (const auto& l2 : leptons) { - if (l1 == l2) continue; - if ((l1.pdgId() + l2.pdgId() == 0) && ((l1.mom() + l2.mom()).mass() < 5.0*GeV)) return false; + for (const Particle& l1 : leptons) { + for (const Particle& l2 : leptons) { + if (isSame(l1, l2)) continue; + if ((l1.pid() + l2.pid() == 0) && ((l1.mom() + l2.mom()).mass() < 5.0*GeV)) return false; } } return true; } - + bool passesLeptonDeltaRRequirements(const Quadruplet& quadruplet) { const auto& leptons = quadruplet.leptonsSortedByPt; - for (const auto& l1 : leptons) { - for (const auto& l2 : leptons) { - if (l1 == l2) continue; + for (const Particle& l1 : leptons) { + for (const Particle& l2 : leptons) { + if (isSame(l1, l2)) continue; // Any lepton flavour: if (deltaR(l1.mom(), l2.mom()) < 0.1) return false; // Different lepton flavour: if ((l1.abspid() - l2.abspid() != 0) && (deltaR(l1.mom(), l2.mom()) < 0.2)) return false; } } return true; } - + Quadruplets formQuadrupletsByChannel(const vector& same_flavour_leptons, vector indices) { Quadruplets quadruplets = {}; for (size_t i = 0; i < indices.size(); ++i) { for (size_t k = i+1; k < indices.size(); ++k) { const auto& pair_i = indices.at(i); const auto& pair_k = indices.at(k); if (indicesOverlap(pair_i, pair_k)) continue; const auto d1 = Dilepton({same_flavour_leptons.at(pair_i.first), same_flavour_leptons.at(pair_i.second)}); const auto d2 = Dilepton({same_flavour_leptons.at(pair_k.first), same_flavour_leptons.at(pair_k.second)}); const auto quadruplet = Quadruplet(d1, d2); if (passesHierarchicalPtRequirements(quadruplet)) quadruplets.push_back(quadruplet); } } return quadruplets; } - - Quadruplets formQuadrupletsByChannel(const vector& electrons, vector e_indices, + + Quadruplets formQuadrupletsByChannel(const vector& electrons, vector e_indices, const vector& muons, vector m_indices) { Quadruplets quadruplets = {}; for (const auto& pair_e : e_indices) { for (const auto& pair_m : m_indices) { const auto d1 = Dilepton({electrons.at(pair_e.first), electrons.at(pair_e.second)}); const auto d2 = Dilepton({muons.at(pair_m.first), muons.at(pair_m.second)}); const auto quadruplet = Quadruplet(d1, d2); if (passesHierarchicalPtRequirements(quadruplet)) quadruplets.push_back(quadruplet); } } return quadruplets; } - - + + Quadruplets getQuadruplets(const vector& electrons, const vector& muons) { const auto oc_electrons_indices = getOppositeChargePairsIndices(electrons); const auto oc_muons_indices = getOppositeChargePairsIndices(muons); const auto electron_quadruplets = formQuadrupletsByChannel(electrons, oc_electrons_indices); const auto muon_quadruplets = formQuadrupletsByChannel(muons, oc_muons_indices); const auto mixed_quadruplets = formQuadrupletsByChannel(electrons, oc_electrons_indices, muons, oc_muons_indices); auto quadruplets = electron_quadruplets; quadruplets.insert(quadruplets.end(), muon_quadruplets.begin(), muon_quadruplets.end()); quadruplets.insert(quadruplets.end(), mixed_quadruplets.begin(), mixed_quadruplets.end()); - + return quadruplets; } - + Quadruplet selectQuadruplet(const Quadruplets& quadruplets) { if (quadruplets.empty()) throw std::logic_error("Expect at least one quadruplet! The user should veto events without quadruplets."); Quadruplets sortedQuadruplets = quadruplets; std::sort(sortedQuadruplets.begin(), sortedQuadruplets.end(), [](const Quadruplet& a, const Quadruplet& b) { return a.distanceFromZMass() < b.distanceFromZMass(); }); return sortedQuadruplets.at(0); } /// Book histograms and initialise projections before the run void init() { const Cut presel = Cuts::abseta < 5 && Cuts::pT > 100*MeV; const FinalState fs(presel); - + // Prompt leptons, photons, neutrinos // Excluding those from tau decay const PromptFinalState photons(presel && Cuts::abspid == PID::PHOTON, false); const PromptFinalState bare_elecs(presel && Cuts::abspid == PID::ELECTRON, false); const PromptFinalState bare_muons(presel && Cuts::abspid == PID::MUON, false); // Baseline lepton and jet declaration const Cut lepton_baseline_cuts = Cuts::abseta < 2.7 && Cuts::pT > 5*GeV; const DressedLeptons elecs = DressedLeptons(photons, bare_elecs, 0.1, lepton_baseline_cuts); const DressedLeptons muons = DressedLeptons(photons, bare_muons, 0.1, lepton_baseline_cuts); declare(elecs, "electrons"); declare(muons, "muons"); - + VetoedFinalState jet_input(fs); jet_input.addVetoOnThisFinalState(elecs); jet_input.addVetoOnThisFinalState(muons); declare(FastJets(jet_input, FastJets::ANTIKT, 0.4), "jets"); // // Book histograms _h["pT_4l"] = bookHisto1D(2, 1, 1); _h["pT_leading_dilepton"] = bookHisto1D(8, 1, 1); _h["pT_subleading_dilepton"] = bookHisto1D(14, 1, 1); _h["pT_lepton1"] = bookHisto1D(20, 1, 1); _h["pT_lepton2"] = bookHisto1D(26, 1, 1); _h["pT_lepton3"] = bookHisto1D(32, 1, 1); _h["pT_lepton4"] = bookHisto1D(38, 1, 1); _h["absy_4l"] = bookHisto1D(44, 1, 1); _h["deltay_dileptons"] = bookHisto1D(50, 1, 1); _h["deltaphi_dileptons"] = bookHisto1D(56, 1, 1); _h["N_jets"] = bookHisto1D(62, 1, 1); _h["N_central_jets"] = bookHisto1D(68, 1, 1); _h["N_jets60"] = bookHisto1D(74, 1, 1); _h["mass_dijet"] = bookHisto1D(80, 1, 1); _h["deltay_dijet"] = bookHisto1D(86, 1, 1); _h["scalarpTsum_jets"] = bookHisto1D(92, 1, 1); _h["abseta_jet1"] = bookHisto1D(98, 1, 1); _h["abseta_jet2"] = bookHisto1D(104, 1, 1); _h["pT_jet1"] = bookHisto1D(110, 1, 1); _h["pT_jet2"] = bookHisto1D(116, 1, 1); } /// Perform the per-event analysis void analyze(Event const & event) { const double weight = event.weight(); - + const auto& baseline_electrons = apply(event, "electrons").dressedLeptons(); const auto& baseline_muons = apply(event, "muons").dressedLeptons(); - + // Form all possible quadruplets passing hierarchical lepton pT cuts const auto quadruplets = getQuadruplets(baseline_electrons, baseline_muons); - + if (quadruplets.empty()) vetoEvent; - + // Select the best quadruplet, the one minimising the total distance from the Z pole mass auto const quadruplet = selectQuadruplet(quadruplets); - + // Event selection on the best quadruplet if (!passesDileptonMinimumMassRequirement(quadruplet)) vetoEvent; if (!passesLeptonDeltaRRequirements(quadruplet)) vetoEvent; if (!inRange(quadruplet.leadingDilepton.momentum().mass(), 66*GeV, 116*GeV)) vetoEvent; if (!inRange(quadruplet.subleadingDilepton.momentum().mass(), 66*GeV, 116*GeV)) vetoEvent; - + // Select jets Jets alljets = apply(event, "jets").jetsByPt(Cuts::pT > 30*GeV); for (const DressedLepton& lep : quadruplet.leptonsSortedByPt) ifilter_discard(alljets, deltaRLess(lep, 0.4)); const Jets jets = alljets; const Jets centralJets = filterBy(jets, Cuts::abseta < 2.4); const Jets pt60Jets = filterBy(jets, Cuts::pT > 60*GeV); - + const auto& leadingDilepton = quadruplet.leadingDilepton.momentum(); const auto& subleadingDilepton = quadruplet.subleadingDilepton.momentum(); - + _h["pT_4l"]->fill((leadingDilepton + subleadingDilepton).pt()/GeV, weight); _h["pT_leading_dilepton"]->fill(leadingDilepton.pt()/GeV, weight); _h["pT_subleading_dilepton"]->fill(subleadingDilepton.pt()/GeV, weight); _h["pT_lepton1"]->fill(quadruplet.leptonsSortedByPt.at(0).pt()/GeV, weight); _h["pT_lepton2"]->fill(quadruplet.leptonsSortedByPt.at(1).pt()/GeV, weight); _h["pT_lepton3"]->fill(quadruplet.leptonsSortedByPt.at(2).pt()/GeV, weight); _h["pT_lepton4"]->fill(quadruplet.leptonsSortedByPt.at(3).pt()/GeV, weight); _h["absy_4l"]->fill((leadingDilepton + subleadingDilepton).absrapidity(), weight); _h["deltay_dileptons"]->fill(fabs(leadingDilepton.rapidity() - subleadingDilepton.rapidity()), weight); _h["deltaphi_dileptons"]->fill(deltaPhi(leadingDilepton, subleadingDilepton)/pi, weight); _h["N_jets"]->fill(jets.size(), weight); _h["N_central_jets"]->fill(centralJets.size(), weight); _h["N_jets60"]->fill(pt60Jets.size(), weight); - + // If at least one jet present if (jets.empty()) vetoEvent; _h["scalarpTsum_jets"]->fill(sum(jets, pT, 0.)/GeV, weight); _h["abseta_jet1"]->fill(jets.front().abseta(), weight); _h["pT_jet1"]->fill(jets.front().pt()/GeV, weight); - + // If at least two jets present - if (jets.size() < 2) vetoEvent; + if (jets.size() < 2) vetoEvent; _h["mass_dijet"]->fill((jets.at(0).mom() + jets.at(1).mom()).mass()/GeV, weight); _h["deltay_dijet"]->fill(fabs(jets.at(0).rapidity() - jets.at(1).rapidity()), weight); _h["abseta_jet2"]->fill(jets.at(1).abseta(), weight); _h["pT_jet2"]->fill(jets.at(1).pt()/GeV, weight); } /// Normalise histograms etc., after the run void finalize() { // Normalise histograms to cross section const double sf = crossSectionPerEvent() / femtobarn; for (map::iterator it = _h.begin(); it != _h.end(); ++it) { scale(it->second, sf); } } //@} - + private: /// @name Histograms //@{ map _h; static constexpr double Z_mass = 91.1876; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2017_I1625109); } diff --git a/analyses/pluginATLAS/ATLAS_2017_I1626105.cc b/analyses/pluginATLAS/ATLAS_2017_I1626105.cc --- a/analyses/pluginATLAS/ATLAS_2017_I1626105.cc +++ b/analyses/pluginATLAS/ATLAS_2017_I1626105.cc @@ -1,130 +1,129 @@ #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 { - /// @brief lepton differential ttbar analysis at 8 TeV + + /// Lepton differential ttbar analysis at 8 TeV class ATLAS_2017_I1626105 : public Analysis { public: - + + DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2017_I1626105); - + + void init() { Cut eta_full = Cuts::abseta < 5.0 && Cuts::pT > 1.0*MeV; // All final state particles const FinalState fs; // Get photons to dress leptons IdentifiedFinalState photons(fs); photons.acceptIdPair(PID::PHOTON); // Projection to find the electrons PromptFinalState prompt_el(Cuts::abspid == PID::ELECTRON, true); DressedLeptons elecs(photons, prompt_el, 0.1, (Cuts::abseta < 2.5) && (Cuts::pT > 25*GeV)); DressedLeptons veto_elecs(photons, prompt_el, 0.1, eta_full, false); declare(elecs, "elecs"); // Projection to find the muons PromptFinalState prompt_mu(Cuts::abspid == PID::MUON, true); DressedLeptons muons(photons, prompt_mu, 0.1, (Cuts::abseta < 2.5) && (Cuts::pT > 25*GeV)); DressedLeptons veto_muons(photons, prompt_mu, 0.1, eta_full, false); declare(muons, "muons"); // Jet clustering. VetoedFinalState vfs; vfs.addVetoOnThisFinalState(veto_elecs); vfs.addVetoOnThisFinalState(veto_muons); FastJets jets(vfs, FastJets::ANTIKT, 0.4); jets.useInvisibles(); declare(jets, "jets"); // Book histograms bookHistos("lep_pt", 1); bookHistos("lep_eta", 3); bookHistos("dilep_pt", 5); bookHistos("dilep_mass", 7); bookHistos("dilep_rap", 9); bookHistos("dilep_dphi", 11); bookHistos("dilep_sumpt", 13); bookHistos("dilep_sumE", 15); } + void analyze(const Event& event) { vector elecs = sortByPt(apply(event, "elecs").dressedLeptons()); - vector muons = sortByPt(apply(event, "muons").dressedLeptons()); + vector muons = sortByPt(apply(event, "muons").dressedLeptons()); Jets jets = apply(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5); - + // Check overlap of jets/leptons. for (const Jet& jet : jets) { - for (const DressedLepton& el : elecs) { - if (deltaR(jet, el) < 0.4) delete el; - } - for (const DressedLepton& mu : muons) { - if (deltaR(jet, mu) < 0.4) delete mu; - } + ifilter_discard(elecs, deltaRLess(jet, 0.4)); + ifilter_discard(muons, deltaRLess(jet, 0.4)); } - if (elecs.empty() || muons.empty()) vetoEvent; - - if (elecs[0].charge() == muons[0].charge()) vetoEvent; - + if (elecs.empty() || muons.empty()) vetoEvent; + if (elecs[0].charge() == muons[0].charge()) vetoEvent; + FourMomentum el = elecs[0].momentum(); FourMomentum mu = muons[0].momentum(); FourMomentum ll = elecs[0].momentum() + muons[0].momentum(); - + // Fill histograms const double weight = event.weight(); fillHistos("lep_pt", el.pT()/GeV, weight); fillHistos("lep_pt", mu.pT()/GeV, weight); fillHistos("lep_eta", el.abseta(), weight); fillHistos("lep_eta", mu.abseta(), weight); fillHistos("dilep_pt", ll.pT()/GeV, weight); fillHistos("dilep_mass", ll.mass()/GeV, weight); fillHistos("dilep_rap", ll.absrap(), weight); fillHistos("dilep_dphi", deltaPhi(el, mu), weight); fillHistos("dilep_sumpt", (el.pT() + mu.pT())/GeV, weight); fillHistos("dilep_sumE", (el.E() + mu.E())/GeV, weight); } void finalize() { // Normalize to cross-section const double sf = crossSection()/femtobarn/sumOfWeights(); for (auto& hist : _h) { const double norm = 1.0 / hist.second->integral(); - // add overflow to last bin + // add overflow to last bin double overflow = hist.second->overflow().effNumEntries(); hist.second->fillBin(hist.second->numBins() - 1, overflow); // histogram normalisation if (hist.first.find("norm") != string::npos) scale(hist.second, norm); else scale(hist.second, sf); } } private: /// @name Histogram helper functions //@{ void bookHistos(const std::string name, unsigned int index) { _h[name] = bookHisto1D(index, 1, 1); _h["norm_" + name] = bookHisto1D(index + 1, 1, 1); } void fillHistos(const std::string name, double value, double weight) { _h[name]->fill(value, weight); _h["norm_" + name]->fill(value, weight); } map _h; //@} }; // Declare the class as a hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2017_I1626105); } diff --git a/analyses/pluginCMS/CMS_2016_I1491950.cc b/analyses/pluginCMS/CMS_2016_I1491950.cc --- a/analyses/pluginCMS/CMS_2016_I1491950.cc +++ b/analyses/pluginCMS/CMS_2016_I1491950.cc @@ -1,505 +1,503 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Tools/ParticleName.hh" #include "Rivet/Tools/ParticleIdUtils.hh" -namespace Rivet -{ - namespace CMS_2016_I1491950_INTERNAL { //< only visible in this compilation unit +namespace Rivet { + + + namespace { //< only visible in this compilation unit /// @brief Special dressed lepton finder /// /// Find dressed leptons by clustering all leptons and photons class SpecialDressedLeptons : public FinalState { public: /// The default constructor. May specify cuts SpecialDressedLeptons(const FinalState& fs, const Cut& cut) : FinalState(cut) { setName("SpecialDressedLeptons"); IdentifiedFinalState ifs(fs); ifs.acceptIdPair(PID::PHOTON); ifs.acceptIdPair(PID::ELECTRON); ifs.acceptIdPair(PID::MUON); addProjection(ifs, "IFS"); addProjection(FastJets(ifs, FastJets::ANTIKT, 0.1), "LeptonJets"); } - + /// Clone on the heap. virtual unique_ptr clone() const { return unique_ptr(new SpecialDressedLeptons(*this)); } - + /// Retrieve the dressed leptons const vector& dressedLeptons() const { return _clusteredLeptons; } - + private: /// Container which stores the clustered lepton objects vector _clusteredLeptons; - + public: void project(const Event& e) { _theParticles.clear(); _clusteredLeptons.clear(); - + vector allClusteredLeptons; - + const Jets jets = applyProjection(e, "LeptonJets").jetsByPt(5.*GeV); foreach (const Jet& jet, jets) { Particle lepCand; for (const Particle& cand : jet.particles()) { const int absPdgId = abs(cand.pdgId()); if (absPdgId == PID::ELECTRON || absPdgId == PID::MUON) { if (cand.pt() > lepCand.pt()) lepCand = cand; } } //Central lepton must be the major component if ((lepCand.pt() < jet.pt()/2.) || (lepCand.pdgId() == 0)) continue; - - DressedLepton lepton = DressedLepton(lepCand); - + + DressedLepton lepton(lepCand); for (const Particle& cand : jet.particles()) { - if (cand == lepCand) continue; - if (cand.pid() != PID::PHOTON) continue; + if (isSame(cand, lepCand)) continue; + if (cand.pid() != PID::PHOTON) continue; lepton.addPhoton(cand, true); } allClusteredLeptons.push_back(lepton); } - + for (const DressedLepton& lepton : allClusteredLeptons) { if (accept(lepton)) { _clusteredLeptons.push_back(lepton); _theParticles.push_back(lepton.constituentLepton()); _theParticles += lepton.constituentPhotons(); } } } }; } - class CMS_2016_I1491950 : public Analysis - { - public: - typedef CMS_2016_I1491950_INTERNAL::SpecialDressedLeptons SpecialDressedLeptons; + class CMS_2016_I1491950 : public Analysis { + public: + /// Constructor CMS_2016_I1491950() : Analysis("CMS_2016_I1491950") { } /// Book histograms and initialise projections before the run void init() { FinalState fs(Cuts::pT > 0. && Cuts::abseta < 6.); PromptFinalState prompt_fs(fs); prompt_fs.acceptMuonDecays(true); prompt_fs.acceptTauDecays(true); - + // Projection for dressed electrons and muons Cut leptonCuts = Cuts::abseta < 2.5 and Cuts::pt > 30.*GeV; SpecialDressedLeptons dressedleptons(prompt_fs, leptonCuts); addProjection(dressedleptons, "DressedLeptons"); - + // Neutrinos IdentifiedFinalState neutrinos(prompt_fs); neutrinos.acceptNeutrinos(); addProjection(neutrinos, "Neutrinos"); - + // Projection for jets VetoedFinalState fsForJets(fs); fsForJets.addVetoOnThisFinalState(dressedleptons); fsForJets.addVetoOnThisFinalState(neutrinos); addProjection(FastJets(fsForJets, FastJets::ANTIKT, 0.4, JetAlg::DECAY_MUONS, JetAlg::DECAY_INVISIBLES), "Jets"); - + //book hists _hist_thadpt = bookHisto1D("d01-x02-y01"); _hist_thady = bookHisto1D("d03-x02-y01"); _hist_tleppt = bookHisto1D("d05-x02-y01"); _hist_tlepy = bookHisto1D("d07-x02-y01"); _hist_ttpt = bookHisto1D("d09-x02-y01"); _hist_tty = bookHisto1D("d13-x02-y01"); _hist_ttm = bookHisto1D("d11-x02-y01"); _hist_njet = bookHisto1D("d15-x02-y01"); _hist_njets_thadpt_1 = bookHisto1D("d17-x02-y01"); _hist_njets_thadpt_2 = bookHisto1D("d18-x02-y01"); _hist_njets_thadpt_3 = bookHisto1D("d19-x02-y01"); _hist_njets_thadpt_4 = bookHisto1D("d20-x02-y01"); _hist_njets_ttpt_1 = bookHisto1D("d22-x02-y01"); _hist_njets_ttpt_2 = bookHisto1D("d23-x02-y01"); _hist_njets_ttpt_3 = bookHisto1D("d24-x02-y01"); _hist_njets_ttpt_4 = bookHisto1D("d25-x02-y01"); _hist_thady_thadpt_1 = bookHisto1D("d27-x02-y01"); _hist_thady_thadpt_2 = bookHisto1D("d28-x02-y01"); _hist_thady_thadpt_3 = bookHisto1D("d29-x02-y01"); _hist_thady_thadpt_4 = bookHisto1D("d30-x02-y01"); _hist_ttm_tty_1 = bookHisto1D("d32-x02-y01"); _hist_ttm_tty_2 = bookHisto1D("d33-x02-y01"); _hist_ttm_tty_3 = bookHisto1D("d34-x02-y01"); _hist_ttm_tty_4 = bookHisto1D("d35-x02-y01"); _hist_ttpt_ttm_1 = bookHisto1D("d37-x02-y01"); _hist_ttpt_ttm_2 = bookHisto1D("d38-x02-y01"); _hist_ttpt_ttm_3 = bookHisto1D("d39-x02-y01"); _hist_ttpt_ttm_4 = bookHisto1D("d40-x02-y01"); _histnorm_thadpt = bookHisto1D("d42-x02-y01"); _histnorm_thady = bookHisto1D("d44-x02-y01"); _histnorm_tleppt = bookHisto1D("d46-x02-y01"); _histnorm_tlepy = bookHisto1D("d48-x02-y01"); _histnorm_ttpt = bookHisto1D("d50-x02-y01"); _histnorm_tty = bookHisto1D("d54-x02-y01"); _histnorm_ttm = bookHisto1D("d52-x02-y01"); _histnorm_njet = bookHisto1D("d56-x02-y01"); _histnorm_njets_thadpt_1 = bookHisto1D("d58-x02-y01"); _histnorm_njets_thadpt_2 = bookHisto1D("d59-x02-y01"); _histnorm_njets_thadpt_3 = bookHisto1D("d60-x02-y01"); _histnorm_njets_thadpt_4 = bookHisto1D("d61-x02-y01"); _histnorm_njets_ttpt_1 = bookHisto1D("d63-x02-y01"); _histnorm_njets_ttpt_2 = bookHisto1D("d64-x02-y01"); _histnorm_njets_ttpt_3 = bookHisto1D("d65-x02-y01"); _histnorm_njets_ttpt_4 = bookHisto1D("d66-x02-y01"); _histnorm_thady_thadpt_1 = bookHisto1D("d68-x02-y01"); _histnorm_thady_thadpt_2 = bookHisto1D("d69-x02-y01"); _histnorm_thady_thadpt_3 = bookHisto1D("d70-x02-y01"); _histnorm_thady_thadpt_4 = bookHisto1D("d71-x02-y01"); _histnorm_ttm_tty_1 = bookHisto1D("d73-x02-y01"); _histnorm_ttm_tty_2 = bookHisto1D("d74-x02-y01"); _histnorm_ttm_tty_3 = bookHisto1D("d75-x02-y01"); _histnorm_ttm_tty_4 = bookHisto1D("d76-x02-y01"); _histnorm_ttpt_ttm_1 = bookHisto1D("d78-x02-y01"); _histnorm_ttpt_ttm_2 = bookHisto1D("d79-x02-y01"); _histnorm_ttpt_ttm_3 = bookHisto1D("d80-x02-y01"); _histnorm_ttpt_ttm_4 = bookHisto1D("d81-x02-y01"); } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = event.weight(); - + // leptons const SpecialDressedLeptons& dressedleptons_proj = applyProjection(event, "DressedLeptons"); std::vector dressedLeptons = dressedleptons_proj.dressedLeptons(); - + if(dressedLeptons.size() != 1) return; - + // neutrinos const Particles neutrinos = applyProjection(event, "Neutrinos").particlesByPt(); _nusum = FourMomentum(0., 0., 0., 0.); for(const Particle& neutrino : neutrinos) { _nusum += neutrino.momentum(); } _wl = _nusum + dressedLeptons[0].momentum(); - + // jets Cut jet_cut = (Cuts::abseta < 2.5) and (Cuts::pT > 25.*GeV); const Jets jets = applyProjection(event, "Jets").jetsByPt(jet_cut); Jets allJets; for (const Jet& jet : jets) { allJets.push_back(jet); - } + } Jets bJets; for (const Jet& jet : allJets) { if (jet.bTagged()) bJets.push_back(jet); } if(bJets.size() < 2 || allJets.size() < 4) return; - + //construct top quark proxies double Kmin = numeric_limits::max(); for(const Jet& itaj : allJets) { for(const Jet& itbj : allJets) { if (itaj.momentum() == itbj.momentum()) continue; FourMomentum wh(itaj.momentum() + itbj.momentum()); for(const Jet& ithbj : bJets) { if(itaj.momentum() == ithbj.momentum() || itbj.momentum() == ithbj.momentum()) continue; FourMomentum th(wh + ithbj.momentum()); for(const Jet& itlbj : bJets) { if(itaj.momentum() == itlbj.momentum() || itbj.momentum() == itlbj.momentum() || ithbj.momentum() == itlbj.momentum()) continue; FourMomentum tl(_wl + itlbj.momentum()); double K = pow(wh.mass() - 80.4, 2) + pow(th.mass() - 172.5, 2) + pow(tl.mass() - 172.5, 2); if(K < Kmin) { Kmin = K; _tl = tl; _th = th; _wh = wh; } } } } } - _hist_thadpt->fill(_th.pt(), weight); + _hist_thadpt->fill(_th.pt(), weight); _hist_thady->fill(abs(_th.rapidity()) , weight); _hist_tleppt->fill(_tl.pt() , weight); _hist_tlepy->fill(abs(_tl.rapidity()) , weight); - _histnorm_thadpt->fill(_th.pt(), weight); + _histnorm_thadpt->fill(_th.pt(), weight); _histnorm_thady->fill(abs(_th.rapidity()) , weight); _histnorm_tleppt->fill(_tl.pt() , weight); _histnorm_tlepy->fill(abs(_tl.rapidity()) , weight); FourMomentum tt(_tl+_th); _hist_ttpt->fill(tt.pt() , weight); _hist_tty->fill(abs(tt.rapidity()) , weight); _hist_ttm->fill(tt.mass() , weight); _hist_njet->fill(min(allJets.size()-4., 4.), weight); _histnorm_ttpt->fill(tt.pt() , weight); _histnorm_tty->fill(abs(tt.rapidity()) , weight); _histnorm_ttm->fill(tt.mass() , weight); _histnorm_njet->fill(min(allJets.size()-4., 4.), weight); if(allJets.size() == 4) { _hist_njets_thadpt_1->fill(_th.pt(), weight); _hist_njets_ttpt_1->fill(tt.pt(), weight); _histnorm_njets_thadpt_1->fill(_th.pt(), weight); _histnorm_njets_ttpt_1->fill(tt.pt(), weight); } else if(allJets.size() == 5) { _hist_njets_thadpt_2->fill(_th.pt(), weight); _hist_njets_ttpt_2->fill(tt.pt(), weight); _histnorm_njets_thadpt_2->fill(_th.pt(), weight); _histnorm_njets_ttpt_2->fill(tt.pt(), weight); } else if(allJets.size() == 6) { _hist_njets_thadpt_3->fill(_th.pt(), weight); _hist_njets_ttpt_3->fill(tt.pt(), weight); _histnorm_njets_thadpt_3->fill(_th.pt(), weight); _histnorm_njets_ttpt_3->fill(tt.pt(), weight); } else //>= 4 jets { _hist_njets_thadpt_4->fill(_th.pt(), weight); _hist_njets_ttpt_4->fill(tt.pt(), weight); _histnorm_njets_thadpt_4->fill(_th.pt(), weight); _histnorm_njets_ttpt_4->fill(tt.pt(), weight); } if(abs(_th.rapidity()) < 0.5) { _hist_thady_thadpt_1->fill(_th.pt(), weight); _histnorm_thady_thadpt_1->fill(_th.pt(), weight); } else if(abs(_th.rapidity()) < 1.0) { _hist_thady_thadpt_2->fill(_th.pt(), weight); _histnorm_thady_thadpt_2->fill(_th.pt(), weight); } else if(abs(_th.rapidity()) < 1.5) { _hist_thady_thadpt_3->fill(_th.pt(), weight); _histnorm_thady_thadpt_3->fill(_th.pt(), weight); } else if(abs(_th.rapidity()) < 2.5) { _hist_thady_thadpt_4->fill(_th.pt(), weight); _histnorm_thady_thadpt_4->fill(_th.pt(), weight); } if(tt.mass() >= 300. && tt.mass() < 450.) { _hist_ttm_tty_1->fill(abs(tt.rapidity()), weight); _histnorm_ttm_tty_1->fill(abs(tt.rapidity()), weight); } else if(tt.mass() >= 450. && tt.mass() < 625.) { _hist_ttm_tty_2->fill(abs(tt.rapidity()), weight); _histnorm_ttm_tty_2->fill(abs(tt.rapidity()), weight); } else if(tt.mass() >= 625. && tt.mass() < 850.) { _hist_ttm_tty_3->fill(abs(tt.rapidity()), weight); _histnorm_ttm_tty_3->fill(abs(tt.rapidity()), weight); } else if(tt.mass() >= 850. && tt.mass() < 2000.) { _hist_ttm_tty_4->fill(abs(tt.rapidity()), weight); _histnorm_ttm_tty_4->fill(abs(tt.rapidity()), weight); } if(tt.pt() < 35.) { _hist_ttpt_ttm_1->fill(tt.mass(), weight); _histnorm_ttpt_ttm_1->fill(tt.mass(), weight); } else if(tt.pt() < 80.) { _hist_ttpt_ttm_2->fill(tt.mass(), weight); _histnorm_ttpt_ttm_2->fill(tt.mass(), weight); } else if(tt.pt() < 140.) { _hist_ttpt_ttm_3->fill(tt.mass(), weight); _histnorm_ttpt_ttm_3->fill(tt.mass(), weight); } else if(tt.pt() < 500.) { _hist_ttpt_ttm_4->fill(tt.mass(), weight); _histnorm_ttpt_ttm_4->fill(tt.mass(), weight); } } /// Normalise histograms etc., after the run void finalize() { scale(_hist_thadpt, crossSection()/sumOfWeights()); scale(_hist_thady, crossSection()/sumOfWeights()); scale(_hist_tleppt, crossSection()/sumOfWeights()); scale(_hist_tlepy, crossSection()/sumOfWeights()); scale(_hist_ttpt, crossSection()/sumOfWeights()); scale(_hist_tty, crossSection()/sumOfWeights()); scale(_hist_ttm, crossSection()/sumOfWeights()); scale(_hist_njet, crossSection()/sumOfWeights()); scale(_hist_njets_thadpt_1, crossSection()/sumOfWeights()); scale(_hist_njets_thadpt_2, crossSection()/sumOfWeights()); scale(_hist_njets_thadpt_3, crossSection()/sumOfWeights()); scale(_hist_njets_thadpt_4, crossSection()/sumOfWeights()); scale(_hist_njets_ttpt_1, crossSection()/sumOfWeights()); scale(_hist_njets_ttpt_2, crossSection()/sumOfWeights()); scale(_hist_njets_ttpt_3, crossSection()/sumOfWeights()); scale(_hist_njets_ttpt_4, crossSection()/sumOfWeights()); scale(_hist_thady_thadpt_1, crossSection()/sumOfWeights()/0.5); scale(_hist_thady_thadpt_2, crossSection()/sumOfWeights()/0.5); scale(_hist_thady_thadpt_3, crossSection()/sumOfWeights()/0.5); scale(_hist_thady_thadpt_4, crossSection()/sumOfWeights()/1.0); scale(_hist_ttm_tty_1, crossSection()/sumOfWeights()/150.); scale(_hist_ttm_tty_2, crossSection()/sumOfWeights()/175.); scale(_hist_ttm_tty_3, crossSection()/sumOfWeights()/225.); scale(_hist_ttm_tty_4, crossSection()/sumOfWeights()/1150.); scale(_hist_ttpt_ttm_1, crossSection()/sumOfWeights()/35.); scale(_hist_ttpt_ttm_2, crossSection()/sumOfWeights()/45.); scale(_hist_ttpt_ttm_3, crossSection()/sumOfWeights()/60.); scale(_hist_ttpt_ttm_4, crossSection()/sumOfWeights()/360.); scale(_histnorm_thadpt, 1./_histnorm_thadpt->sumW(false)); scale(_histnorm_thady, 1./_histnorm_thady->sumW(false)); scale(_histnorm_tleppt, 1./_histnorm_tleppt->sumW(false)); scale(_histnorm_tlepy, 1./_histnorm_tlepy->sumW(false)); scale(_histnorm_ttpt, 1./_histnorm_ttpt->sumW(false)); scale(_histnorm_tty, 1./_histnorm_tty->sumW(false)); scale(_histnorm_ttm, 1./_histnorm_ttm->sumW(false)); scale(_histnorm_njet, 1./_histnorm_njet->sumW(false)); double sum_njets_thadpt = _histnorm_njets_thadpt_1->sumW(false) + _histnorm_njets_thadpt_2->sumW(false) + _histnorm_njets_thadpt_3->sumW(false) + _histnorm_njets_thadpt_4->sumW(false); scale(_histnorm_njets_thadpt_1, 1./sum_njets_thadpt); scale(_histnorm_njets_thadpt_2, 1./sum_njets_thadpt); scale(_histnorm_njets_thadpt_3, 1./sum_njets_thadpt); scale(_histnorm_njets_thadpt_4, 1./sum_njets_thadpt); double sum_njets_ttpt = _histnorm_njets_ttpt_1->sumW(false) + _histnorm_njets_ttpt_2->sumW(false) + _histnorm_njets_ttpt_3->sumW(false) + _histnorm_njets_ttpt_4->sumW(false); scale(_histnorm_njets_ttpt_1, 1./sum_njets_ttpt); scale(_histnorm_njets_ttpt_2, 1./sum_njets_ttpt); scale(_histnorm_njets_ttpt_3, 1./sum_njets_ttpt); scale(_histnorm_njets_ttpt_4, 1./sum_njets_ttpt); double sum_thady_thadpt = _histnorm_thady_thadpt_1->sumW(false) + _histnorm_thady_thadpt_2->sumW(false) + _histnorm_thady_thadpt_3->sumW(false) + _histnorm_thady_thadpt_4->sumW(false); scale(_histnorm_thady_thadpt_1, 1./sum_thady_thadpt/0.5); scale(_histnorm_thady_thadpt_2, 1./sum_thady_thadpt/0.5); scale(_histnorm_thady_thadpt_3, 1./sum_thady_thadpt/0.5); scale(_histnorm_thady_thadpt_4, 1./sum_thady_thadpt/1.0); double sum_ttm_tty = _histnorm_ttm_tty_1->sumW(false) + _histnorm_ttm_tty_2->sumW(false) + _histnorm_ttm_tty_3->sumW(false) + _histnorm_ttm_tty_4->sumW(false); scale(_histnorm_ttm_tty_1, 1./sum_ttm_tty/150.); scale(_histnorm_ttm_tty_2, 1./sum_ttm_tty/175.); scale(_histnorm_ttm_tty_3, 1./sum_ttm_tty/225.); scale(_histnorm_ttm_tty_4, 1./sum_ttm_tty/1150.); double sum_ttpt_ttm = _histnorm_ttpt_ttm_1->sumW(false) + _histnorm_ttpt_ttm_2->sumW(false) + _histnorm_ttpt_ttm_3->sumW(false) + _histnorm_ttpt_ttm_4->sumW(false); scale(_histnorm_ttpt_ttm_1, 1./sum_ttpt_ttm/35.); scale(_histnorm_ttpt_ttm_2, 1./sum_ttpt_ttm/45.); scale(_histnorm_ttpt_ttm_3, 1./sum_ttpt_ttm/60.); scale(_histnorm_ttpt_ttm_4, 1./sum_ttpt_ttm/360.); } private: FourMomentum _tl; FourMomentum _th; FourMomentum _wl; FourMomentum _wh; FourMomentum _nusum; Histo1DPtr _hist_thadpt; Histo1DPtr _hist_thady; Histo1DPtr _hist_tleppt; Histo1DPtr _hist_tlepy; Histo1DPtr _hist_ttpt; Histo1DPtr _hist_tty; Histo1DPtr _hist_ttm; Histo1DPtr _hist_njet; Histo1DPtr _hist_njets_thadpt_1; Histo1DPtr _hist_njets_thadpt_2; Histo1DPtr _hist_njets_thadpt_3; Histo1DPtr _hist_njets_thadpt_4; Histo1DPtr _hist_njets_ttpt_1; Histo1DPtr _hist_njets_ttpt_2; Histo1DPtr _hist_njets_ttpt_3; Histo1DPtr _hist_njets_ttpt_4; Histo1DPtr _hist_thady_thadpt_1; Histo1DPtr _hist_thady_thadpt_2; Histo1DPtr _hist_thady_thadpt_3; Histo1DPtr _hist_thady_thadpt_4; Histo1DPtr _hist_ttm_tty_1; Histo1DPtr _hist_ttm_tty_2; Histo1DPtr _hist_ttm_tty_3; Histo1DPtr _hist_ttm_tty_4; Histo1DPtr _hist_ttpt_ttm_1; Histo1DPtr _hist_ttpt_ttm_2; Histo1DPtr _hist_ttpt_ttm_3; Histo1DPtr _hist_ttpt_ttm_4; Histo1DPtr _histnorm_thadpt; Histo1DPtr _histnorm_thady; Histo1DPtr _histnorm_tleppt; Histo1DPtr _histnorm_tlepy; Histo1DPtr _histnorm_ttpt; Histo1DPtr _histnorm_tty; Histo1DPtr _histnorm_ttm; Histo1DPtr _histnorm_njet; Histo1DPtr _histnorm_njets_thadpt_1; Histo1DPtr _histnorm_njets_thadpt_2; Histo1DPtr _histnorm_njets_thadpt_3; Histo1DPtr _histnorm_njets_thadpt_4; Histo1DPtr _histnorm_njets_ttpt_1; Histo1DPtr _histnorm_njets_ttpt_2; Histo1DPtr _histnorm_njets_ttpt_3; Histo1DPtr _histnorm_njets_ttpt_4; Histo1DPtr _histnorm_thady_thadpt_1; Histo1DPtr _histnorm_thady_thadpt_2; Histo1DPtr _histnorm_thady_thadpt_3; Histo1DPtr _histnorm_thady_thadpt_4; Histo1DPtr _histnorm_ttm_tty_1; Histo1DPtr _histnorm_ttm_tty_2; Histo1DPtr _histnorm_ttm_tty_3; Histo1DPtr _histnorm_ttm_tty_4; Histo1DPtr _histnorm_ttpt_ttm_1; Histo1DPtr _histnorm_ttpt_ttm_2; Histo1DPtr _histnorm_ttpt_ttm_3; - Histo1DPtr _histnorm_ttpt_ttm_4; - + Histo1DPtr _histnorm_ttpt_ttm_4; + }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_2016_I1491950); } - diff --git a/analyses/pluginCMS/CMS_2016_PAS_TOP_15_006.cc b/analyses/pluginCMS/CMS_2016_PAS_TOP_15_006.cc --- a/analyses/pluginCMS/CMS_2016_PAS_TOP_15_006.cc +++ b/analyses/pluginCMS/CMS_2016_PAS_TOP_15_006.cc @@ -1,186 +1,182 @@ #include "Rivet/Analysis.hh" -#include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" -#include "Rivet/Tools/ParticleName.hh" -#include "Rivet/Tools/ParticleIdUtils.hh" namespace Rivet { - namespace CMS_2016_PAS_TOP_15_006_INTERNAL { //< only visible in this compilation unit + namespace { //< only visible in this compilation unit /// @brief Special dressed lepton finder /// /// Find dressed leptons by clustering all leptons and photons class SpecialDressedLeptons : public FinalState { public: /// Constructor SpecialDressedLeptons(const FinalState& fs, const Cut& cut) : FinalState(cut) { setName("SpecialDressedLeptons"); IdentifiedFinalState ifs(fs); ifs.acceptIdPair(PID::PHOTON); ifs.acceptIdPair(PID::ELECTRON); ifs.acceptIdPair(PID::MUON); addProjection(ifs, "IFS"); addProjection(FastJets(ifs, FastJets::ANTIKT, 0.1), "LeptonJets"); } /// Clone on the heap virtual unique_ptr clone() const { return unique_ptr(new SpecialDressedLeptons(*this)); } /// Retrieve the dressed leptons const vector& dressedLeptons() const { return _clusteredLeptons; } /// Perform the calculation void project(const Event& e) { _theParticles.clear(); _clusteredLeptons.clear(); vector allClusteredLeptons; const Jets jets = applyProjection(e, "LeptonJets").jetsByPt(5*GeV); for (const Jet& jet : jets) { Particle lepCand; for (const Particle& cand : jet.particles()) { const int absPdgId = cand.abspid(); if (absPdgId == PID::ELECTRON || absPdgId == PID::MUON) { if (cand.pt() > lepCand.pt()) lepCand = cand; } } // Central lepton must be the major component if ((lepCand.pt() < jet.pt()/2.) || (lepCand.pdgId() == 0)) continue; DressedLepton lepton = DressedLepton(lepCand); for (const Particle& cand : jet.particles()) { - if (cand == lepCand) continue; + if (isSame(cand, lepCand)) continue; lepton.addPhoton(cand, true); } allClusteredLeptons.push_back(lepton); } for (const DressedLepton& lepton : allClusteredLeptons) { if (accept(lepton)) { _clusteredLeptons.push_back(lepton); _theParticles.push_back(lepton.constituentLepton()); _theParticles += lepton.constituentPhotons(); } } } + private: /// Container which stores the clustered lepton objects vector _clusteredLeptons; }; } /// Jet multiplicity in lepton+jets ttbar at 8 TeV class CMS_2016_PAS_TOP_15_006 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(CMS_2016_PAS_TOP_15_006); - typedef CMS_2016_PAS_TOP_15_006_INTERNAL::SpecialDressedLeptons SpecialDressedLeptons; - /// @name Analysis methods //@{ /// Set up projections and book histograms void init() { // Complete final state FinalState fs; Cut superLooseLeptonCuts = Cuts::pt > 5*GeV; SpecialDressedLeptons dressedleptons(fs, superLooseLeptonCuts); addProjection(dressedleptons, "DressedLeptons"); // Projection for jets VetoedFinalState fsForJets(fs); fsForJets.addVetoOnThisFinalState(dressedleptons); addProjection(FastJets(fsForJets, FastJets::ANTIKT, 0.5), "Jets"); // Booking of histograms _normedElectronMuonHisto = bookHisto1D("normedElectronMuonHisto", 7, 3.5, 10.5, "Normalized differential cross section in lepton+jets channel", "Jet multiplicity", "Normed units"); _absXSElectronMuonHisto = bookHisto1D("absXSElectronMuonHisto", 7, 3.5, 10.5, "Differential cross section in lepton+jets channel", "Jet multiplicity", "pb"); } /// Per-event analysis void analyze(const Event& event) { const double weight = event.weight(); // Select ttbar -> lepton+jets const SpecialDressedLeptons& dressedleptons = applyProjection(event, "DressedLeptons"); vector selleptons; for (const DressedLepton& dressedlepton : dressedleptons.dressedLeptons()) { // Select good leptons if (dressedlepton.pT() > 30*GeV && dressedlepton.abseta() < 2.4) selleptons += dressedlepton.mom(); // Veto loose leptons else if (dressedlepton.pT() > 15*GeV && dressedlepton.abseta() < 2.5) vetoEvent; } if (selleptons.size() != 1) vetoEvent; // Identify hardest tight lepton const FourMomentum lepton = selleptons[0]; // Jets const FastJets& jets = applyProjection(event, "Jets"); const Jets jets30 = jets.jetsByPt(30*GeV); int nJets = 0, nBJets = 0; for (const Jet& jet : jets30) { if (jet.abseta() > 2.5) continue; if (deltaR(jet.momentum(), lepton) < 0.5) continue; nJets += 1; if (jet.bTagged(Cuts::pT > 5*GeV)) nBJets += 1; } // Require >= 4 resolved jets, of which two must be b-tagged if (nJets < 4 || nBJets < 2) vetoEvent; // Fill histograms _normedElectronMuonHisto->fill(min(nJets, 10), weight); _absXSElectronMuonHisto ->fill(min(nJets, 10), weight); } void finalize() { const double ttbarXS = !std::isnan(crossSectionPerEvent()) ? crossSection() : 252.89*picobarn; if (std::isnan(crossSectionPerEvent())) MSG_INFO("No valid cross-section given, using NNLO (arXiv:1303.6254; sqrt(s)=8 TeV, m_t=172.5 GeV): " << ttbarXS/picobarn << " pb"); const double xsPerWeight = ttbarXS/picobarn / sumOfWeights(); scale(_absXSElectronMuonHisto, xsPerWeight); normalize(_normedElectronMuonHisto); } //@} private: /// Histograms Histo1DPtr _normedElectronMuonHisto, _absXSElectronMuonHisto; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_2016_PAS_TOP_15_006); }