Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F7879768
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
41 KB
Subscribers
None
View Options
diff --git a/src/Analyses/ATLAS_2010_S8914702.cc b/src/Analyses/ATLAS_2010_S8914702.cc
--- a/src/Analyses/ATLAS_2010_S8914702.cc
+++ b/src/Analyses/ATLAS_2010_S8914702.cc
@@ -1,202 +1,202 @@
// -*- C++ -*-
#include <iostream>
#include <sstream>
#include <string>
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/LeadingParticlesFinalState.hh"
#include "Rivet/Jet.hh"
#include "Rivet/Projections/FastJets.hh"
#include "fastjet/internal/base.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/AreaDefinition.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include "fastjet/PseudoJet.hh"
namespace Rivet {
class ATLAS_2010_S8914702 : public Analysis {
public:
/// Constructor
ATLAS_2010_S8914702()
: Analysis("ATLAS_2010_S8914702")
{
_eta_bins.push_back( 0.00);
_eta_bins.push_back( 0.60);
_eta_bins.push_back( 1.37);
_eta_bins.push_back( 1.52);
_eta_bins.push_back( 1.81);
_eta_bins_areaoffset.push_back(0.0);
_eta_bins_areaoffset.push_back(1.5);
_eta_bins_areaoffset.push_back(3.0);
}
/// Book histograms and initialise projections before the run
void init() {
FinalState fs;
addProjection(fs, "FS");
FastJets fj(fs, FastJets::KT, 0.5);
_area_def = new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec());
fj.useJetArea(_area_def);
addProjection(fj, "KtJetsD05");
LeadingParticlesFinalState photonfs(FinalState(-1.81, 1.81, 15.0*GeV));
photonfs.addParticleId(PID::PHOTON);
addProjection(photonfs, "LeadingPhoton");
int hist_bin = 0;
for (int i = 0; i < (int)_eta_bins.size()-1; ++i) {
if (fabs(_eta_bins[i] - 1.37) < .0001) continue;
_h_Et_photon[i] = bookHisto1D(1, 1, hist_bin+1);
hist_bin += 1;
}
}
int getEtaBin(double eta_w, bool area_eta) const {
double eta = fabs(eta_w);
int v_iter = 0;
if (!area_eta) {
for (v_iter=0; v_iter < (int)_eta_bins.size()-1; ++v_iter) {
if (eta >= _eta_bins.at(v_iter) && eta < _eta_bins.at(v_iter+1)) break;
}
return min(v_iter,(int)_eta_bins.size()-2);
} else {
for (v_iter=0; v_iter < (int)_eta_bins_areaoffset.size()-1; ++v_iter) {
if (eta >= _eta_bins_areaoffset.at(v_iter) && eta < _eta_bins_areaoffset.at(v_iter+1)) break;
}
return v_iter;
}
}
/// Perform the per-event analysis
void analyze(const Event& event) {
const double weight = event.weight();
Particles photons = applyProjection<LeadingParticlesFinalState>(event, "LeadingPhoton").particles();
if (photons.size() != 1) {
vetoEvent;
}
FourMomentum leadingPhoton = photons[0].momentum();
double eta_P = leadingPhoton.eta();
double phi_P = leadingPhoton.phi();
if(fabs(eta_P)>=1.37 && fabs(eta_P)<1.52){
vetoEvent;
}
int eta_bin = getEtaBin(eta_P,false);
Particles fs = applyProjection<FinalState>(event, "FS").particles();
FourMomentum mom_in_EtCone;
foreach (const Particle& p, fs) {
// check if it's in the cone of .4
if (deltaR(eta_P, phi_P, p.eta(), p.phi()) >= 0.4) continue;
// check if it's in the 5x7 central core
if (fabs(eta_P-p.eta()) < .025*7.0*0.5 &&
- fabs(phi_P-p.phi()) < (PI/128.)*5.0*0.5) continue;
+ fabs(phi_P-p.phi()) < (PI/128.)*7.0*0.5) continue;
mom_in_EtCone += p.momentum();
}
MSG_DEBUG("Done with initial EtCone.");
// Now compute the median energy density
_ptDensity.clear();
_sigma.clear();
_Njets.clear();
std::vector< std::vector<double> > ptDensities;
std::vector<double> emptyVec;
ptDensities.assign(_eta_bins_areaoffset.size()-1,emptyVec);
FastJets fast_jets = applyProjection<FastJets>(event, "KtJetsD05");
const fastjet::ClusterSequenceArea* clust_seq_area = fast_jets.clusterSeqArea();
foreach (const fastjet::PseudoJet& jet, fast_jets.pseudoJets(0.0*GeV)) {
//const double y = jet.absrap();
const double eta = fabs(jet.eta());
const double pt = fabs(jet.perp());
// Get the cluster sequence
double area = clust_seq_area->area(jet);
if (area > 10e-4 && fabs(eta)<_eta_bins_areaoffset[_eta_bins_areaoffset.size()-1]) {
ptDensities.at(getEtaBin(fabs(eta),true)).push_back(pt/area);
}
}
for (int b = 0; b < (int)_eta_bins_areaoffset.size()-1; ++b) {
double median = 0.0;
double sigma = 0.0;
int Njets = 0;
if (ptDensities[b].size() > 0) {
std::sort(ptDensities[b].begin(), ptDensities[b].end());
int nDens = ptDensities[b].size();
if (nDens % 2 == 0) {
median = (ptDensities[b][nDens/2]+ptDensities[b][(nDens-2)/2])/2;
} else {
median = ptDensities[b][(nDens-1)/2];
}
sigma = ptDensities[b][(int)(.15865*nDens)];
Njets = nDens;
}
_ptDensity.push_back(median);
_sigma.push_back(sigma);
_Njets.push_back(Njets);
}
// Now figure out the correction
float EtCone_area = PI*.4*.4 - (7.0*.025)*(5.0*PI/128.);
float correction = _ptDensity[getEtaBin(eta_P,true)]*EtCone_area;
MSG_DEBUG("Jet area correction done.");
// Shouldn't need to subtract photon
// NB. Using expected cut at hadron/particle level, not cut at reco level
if (mom_in_EtCone.Et() - correction/*-leadingPhoton.Et()*/ > 4.0*GeV) {
vetoEvent;
}
MSG_DEBUG("Passed isolation cut.");
_h_Et_photon[eta_bin]->fill(leadingPhoton.Et(), weight);
}
/// Normalise histograms etc., after the run
void finalize() {
for (int i = 0; i < (int)_eta_bins.size()-1; ++i) {
if (fabs(_eta_bins[i] - 1.37) < .0001) continue;
scale(_h_Et_photon[i], crossSection()/sumOfWeights());
}
}
private:
Histo1DPtr _h_Et_photon[6];
fastjet::AreaDefinition* _area_def;
std::vector<float> _eta_bins;
std::vector<float> _eta_bins_areaoffset;
std::vector<float> _ptDensity;
std::vector<float> _sigma;
std::vector<float> _Njets;
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(ATLAS_2010_S8914702);
}
diff --git a/src/Analyses/ATLAS_2011_I921594.cc b/src/Analyses/ATLAS_2011_I921594.cc
--- a/src/Analyses/ATLAS_2011_I921594.cc
+++ b/src/Analyses/ATLAS_2011_I921594.cc
@@ -1,145 +1,145 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/LeadingParticlesFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
namespace Rivet {
/// Inclusive isolated prompt photon analysis with full 2010 LHC data
class ATLAS_2011_I921594 : public Analysis {
public:
/// Constructor
ATLAS_2011_I921594()
: Analysis("ATLAS_2011_I921594")
{
_eta_bins.push_back( 0.00);
_eta_bins.push_back( 0.60);
_eta_bins.push_back( 1.37);
_eta_bins.push_back( 1.52);
_eta_bins.push_back( 1.81);
_eta_bins.push_back( 2.37);
_eta_bins_areaoffset.push_back(0.0);
_eta_bins_areaoffset.push_back(1.5);
_eta_bins_areaoffset.push_back(3.0);
}
/// Book histograms and initialise projections before the run
void init() {
FinalState fs;
addProjection(fs, "FS");
// Consider the final state jets for the energy density calculation
FastJets fj(fs, FastJets::KT, 0.5);
/// @todo Memory leak! (Once per run, not serious)
_area_def = new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec()); ///< @todo FastJets should deep copy the pointer if possible
fj.useJetArea(_area_def);
addProjection(fj, "KtJetsD05");
// Consider the leading pt photon with |eta|<2.37 and pT>45 GeV
LeadingParticlesFinalState photonfs(FinalState(-2.37, 2.37, 45*GeV));
photonfs.addParticleId(PID::PHOTON);
addProjection(photonfs, "LeadingPhoton");
// Book the dsigma/dEt (in eta bins) histograms
for (size_t i = 0; i < _eta_bins.size()-1; i++) {
if (fuzzyEquals(_eta_bins[i], 1.37)) continue; // skip this bin
_h_Et_photon[i] = bookHisto1D(1, 1, i+1);
}
}
/// Return eta bin for either dsigma/dET histogram (area_eta=false) or energy density correction (area_eta=true)
size_t _getEtaBin(double eta_w, bool area_eta) const {
const double eta = fabs(eta_w);
if (!area_eta) {
return binIndex(eta, _eta_bins);
} else {
return binIndex(eta, _eta_bins_areaoffset);
}
}
/// Perform the per-event analysis
void analyze(const Event& event) {
// Retrieve leading photon
const Particles& photons = applyProjection<LeadingParticlesFinalState>(event, "LeadingPhoton").particles();
if (photons.size() != 1) vetoEvent;
const Particle& leadingPhoton = photons[0];
// Veto events with photon in ECAL crack
if (inRange(leadingPhoton.abseta(), 1.37, 1.52)) vetoEvent;
// Compute isolation energy in cone of radius .4 around photon (all particles)
FourMomentum mom_in_EtCone;
Particles fs = applyProjection<FinalState>(event, "FS").particles();
foreach (const Particle& p, fs) {
// Check if it's outside the cone of 0.4
if (deltaR(leadingPhoton, p) >= 0.4) continue;
// Don't count particles in the 5x7 central core
if (deltaEta(leadingPhoton, p) < .025*7.0*0.5 &&
- deltaPhi(leadingPhoton, p) < (PI/128.)*5.0*0.5) continue;
+ deltaPhi(leadingPhoton, p) < (PI/128.)*7.0*0.5) continue;
// Increment isolation energy
mom_in_EtCone += p.momentum();
}
// Get the area-filtered jet inputs for computing median energy density, etc.
vector<double> ptDensity, ptSigma, nJets;
vector< vector<double> > ptDensities(_eta_bins_areaoffset.size()-1);
FastJets fast_jets = applyProjection<FastJets>(event, "KtJetsD05");
const fastjet::ClusterSequenceArea* clust_seq_area = fast_jets.clusterSeqArea();
foreach (const Jet& jet, fast_jets.jets()) {
const double area = clust_seq_area->area(jet);
if (area > 10e-4 && jet.abseta() < _eta_bins_areaoffset.back())
ptDensities.at( _getEtaBin(jet.abseta(), true) ).push_back(jet.pT()/area);
}
// Compute the median energy density, etc.
for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; b++) {
const int njets = ptDensities[b].size();
const double ptmedian = (njets > 0) ? median(ptDensities[b]) : 0;
const double ptsigma = (njets > 0) ? ptDensities[b][(size_t)(0.15865*njets)] : 0;
nJets.push_back(njets);
ptDensity.push_back(ptmedian);
ptSigma.push_back(ptsigma);
}
// Compute the isolation energy correction (cone area*energy density)
const double etCone_area = PI*sqr(0.4) - (7.0*.025)*(5.0*PI/128.);
const double correction = ptDensity[_getEtaBin(leadingPhoton.abseta(), true)]*etCone_area;
// Apply isolation cut on area-corrected value
if (mom_in_EtCone.Et() - correction > 4*GeV) vetoEvent;
// Fill histograms
const size_t eta_bin = _getEtaBin(leadingPhoton.abseta(), false);
_h_Et_photon[eta_bin]->fill(leadingPhoton.Et(), event.weight());
}
/// Normalise histograms etc., after the run
void finalize() {
for (size_t i = 0; i < _eta_bins.size()-1; i++) {
if (fuzzyEquals(_eta_bins[i], 1.37)) continue;
scale(_h_Et_photon[i], crossSection()/picobarn/sumOfWeights());
}
}
private:
Histo1DPtr _h_Et_photon[5];
fastjet::AreaDefinition* _area_def;
vector<double> _eta_bins, _eta_bins_areaoffset;
};
DECLARE_RIVET_PLUGIN(ATLAS_2011_I921594);
}
diff --git a/src/Analyses/ATLAS_2011_S9120807.cc b/src/Analyses/ATLAS_2011_S9120807.cc
--- a/src/Analyses/ATLAS_2011_S9120807.cc
+++ b/src/Analyses/ATLAS_2011_S9120807.cc
@@ -1,223 +1,223 @@
// -*- C++ -*-
#include <iostream>
#include <sstream>
#include <string>
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/Jet.hh"
#include "Rivet/Projections/FastJets.hh"
#include "fastjet/internal/base.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/AreaDefinition.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include "fastjet/PseudoJet.hh"
namespace Rivet {
/// @brief Measurement of isolated diphoton + X differential cross-sections
///
/// Inclusive isolated gamma gamma cross-sections, differential in M(gg), pT(gg),
/// dphi(gg)
///
/// @author Giovanni Marchiori
class ATLAS_2011_S9120807 : public Analysis {
public:
/// Constructor
ATLAS_2011_S9120807()
: Analysis("ATLAS_2011_S9120807")
{
_eta_bins_areaoffset.push_back(0.0);
_eta_bins_areaoffset.push_back(1.5);
_eta_bins_areaoffset.push_back(3.0);
}
public:
/// Book histograms and initialise projections before the run
void init() {
FinalState fs;
addProjection(fs, "FS");
FastJets fj(fs, FastJets::KT, 0.5);
_area_def = new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec());
fj.useJetArea(_area_def);
addProjection(fj, "KtJetsD05");
IdentifiedFinalState photonfs(Cuts::abseta < 2.37 && Cuts::pT > 16*GeV);
photonfs.acceptId(PID::PHOTON);
addProjection(photonfs, "Photon");
_h_M = bookHisto1D(1, 1, 1);
_h_pT = bookHisto1D(2, 1, 1);
_h_dPhi = bookHisto1D(3, 1, 1);
}
/// @todo Prefer to use Rivet::binIndex()
size_t getEtaBin(double eta_w) const {
const double aeta = fabs(eta_w);
size_t v_iter = 0;
for (; v_iter+1 < _eta_bins_areaoffset.size(); ++v_iter) {
if (inRange(aeta, _eta_bins_areaoffset[v_iter], _eta_bins_areaoffset[v_iter+1])) break;
}
return v_iter;
}
/// Perform the per-event analysis
void analyze(const Event& event) {
const double weight = event.weight();
///
/// require at least 2 photons in final state
///
Particles photons = applyProjection<IdentifiedFinalState>(event, "Photon").particlesByPt();
if (photons.size() < 2) {
vetoEvent;
}
///
/// compute the median energy density
///
std::vector<double> _ptDensity;
std::vector< std::vector<double> > ptDensities;
std::vector<double> emptyVec;
ptDensities.assign(_eta_bins_areaoffset.size()-1, emptyVec);
const fastjet::ClusterSequenceArea* clust_seq_area = applyProjection<FastJets>(event, "KtJetsD05").clusterSeqArea();
foreach (const fastjet::PseudoJet& jet, applyProjection<FastJets>(event, "KtJetsD05").pseudoJets(0.0*GeV)) {
double eta = fabs(jet.eta());
double pt = fabs(jet.perp());
/// get the cluster sequence
double area = clust_seq_area->area(jet);
if(area > 10e-4 && fabs(eta) < _eta_bins_areaoffset[_eta_bins_areaoffset.size()-1]) {
ptDensities.at(getEtaBin(fabs(eta))).push_back(pt/area);
}
}
for(int b=0; b<(int)_eta_bins_areaoffset.size()-1; b++) {
double median = 0.0;
if(ptDensities[b].size() > 0) {
std::sort(ptDensities[b].begin(), ptDensities[b].end());
int nDens = ptDensities[b].size();
if( nDens%2 == 0 )
median = (ptDensities[b][nDens/2] + ptDensities[b][(nDens-2)/2]) / 2;
else
median = ptDensities[b][(nDens-1)/2];
}
_ptDensity.push_back(median);
}
///
/// Loop over photons and fill vector of isolated ones
///
Particles isolated_photons;
foreach (const Particle& photon, photons) {
///
/// remove photons in crack
///
double eta_P = photon.eta();
if (fabs(eta_P)>=1.37 && fabs(eta_P)<1.52) continue;
double phi_P = photon.phi();
///
/// compute isolation
///
/// std EtCone
Particles fs = applyProjection<FinalState>(event, "FS").particles();
FourMomentum mom_in_EtCone;
foreach (const Particle& p, fs) {
/// check if it's in the cone of .4
if (deltaR(eta_P, phi_P, p.eta(), p.phi()) >= 0.4) continue;
/// check if it's in the 5x7 central core
if (fabs(eta_P-p.eta()) < .025*7.0*0.5 &&
- fabs(phi_P-p.phi()) < (M_PI/128.)*5.0*0.5) continue;
+ fabs(phi_P-p.phi()) < (M_PI/128.)*7.0*0.5) continue;
mom_in_EtCone += p.momentum();
}
/// now figure out the correction (area*density)
double EtCone_area = M_PI*.4*.4 - (7.0*.025)*(5.0*M_PI/128.);
double correction = _ptDensity[getEtaBin(eta_P)]*EtCone_area;
/// shouldn't need to subtract photon
/// note: using expected cut at hadron/particle level, not cut at reco level
if (mom_in_EtCone.Et()-correction > 4.0*GeV) {
continue;
}
/// add photon to list of isolated ones
isolated_photons.push_back(photon);
}
///
/// require at least two isolated photons
///
if (isolated_photons.size() < 2) {
vetoEvent;
}
///
/// select leading pT pair
///
std::sort(isolated_photons.begin(), isolated_photons.end(), cmpMomByPt);
FourMomentum y1=isolated_photons[0].momentum();
FourMomentum y2=isolated_photons[1].momentum();
///
/// require the two photons to be separated (dR>0.4)
///
if (deltaR(y1, y2)<0.4) {
vetoEvent;
}
FourMomentum yy = y1+y2;
double Myy = yy.mass()/GeV;
double pTyy = yy.pT()/GeV;
double dPhiyy = deltaPhi(y1.phi(), y2.phi());
_h_M->fill(Myy, weight);
_h_pT->fill(pTyy, weight);
_h_dPhi->fill(dPhiyy, weight);
}
/// Normalise histograms etc., after the run
void finalize() {
scale(_h_M, crossSection()/sumOfWeights());
scale(_h_pT, crossSection()/sumOfWeights());
scale(_h_dPhi, crossSection()/sumOfWeights());
}
private:
Histo1DPtr _h_M;
Histo1DPtr _h_pT;
Histo1DPtr _h_dPhi;
fastjet::AreaDefinition* _area_def;
std::vector<double> _eta_bins_areaoffset;
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(ATLAS_2011_S9120807);
}
diff --git a/src/Analyses/ATLAS_2012_I1093738.cc b/src/Analyses/ATLAS_2012_I1093738.cc
--- a/src/Analyses/ATLAS_2012_I1093738.cc
+++ b/src/Analyses/ATLAS_2012_I1093738.cc
@@ -1,288 +1,288 @@
// -*- C++ -*-
#include <iostream>
#include <sstream>
#include <string>
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/LeadingParticlesFinalState.hh"
#include "Rivet/Projections/VetoedFinalState.hh"
#include "Rivet/Jet.hh"
#include "Rivet/Projections/FastJets.hh"
#include "fastjet/internal/base.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/AreaDefinition.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include "fastjet/PseudoJet.hh"
namespace Rivet {
/// @brief Measurement of isolated gamma + jet + X differential cross-sections
///
/// Inclusive isolated gamma + jet cross-sections, differential in pT(gamma), for
/// various photon and jet rapidity configurations.
///
/// @author Giovanni Marchiori
class ATLAS_2012_I1093738 : public Analysis {
public:
// Constructor
ATLAS_2012_I1093738()
: Analysis("ATLAS_2012_I1093738")
{
_eta_bins_ph.push_back(0.0);
_eta_bins_ph.push_back(1.37);
_eta_bins_ph.push_back(1.52);
_eta_bins_ph.push_back(2.37);
_eta_bins_jet.push_back(0.0);
_eta_bins_jet.push_back(1.2);
_eta_bins_jet.push_back(2.8);
_eta_bins_jet.push_back(4.4);
_eta_bins_areaoffset.push_back(0.0);
_eta_bins_areaoffset.push_back(1.5);
_eta_bins_areaoffset.push_back(3.0);
}
public:
// Book histograms and initialise projections before the run
void init() {
// Final state
FinalState fs;
addProjection(fs, "FS");
// Voronoi eta-phi tassellation with KT jets, for ambient energy density calculation
FastJets fj(fs, FastJets::KT, 0.5);
_area_def = new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec());
fj.useJetArea(_area_def);
addProjection(fj, "KtJetsD05");
// Leading photon
LeadingParticlesFinalState photonfs(FinalState(-1.37, 1.37, 25.0*GeV));
photonfs.addParticleId(PID::PHOTON);
addProjection(photonfs, "LeadingPhoton");
// FS excluding the leading photon
VetoedFinalState vfs(fs);
vfs.addVetoOnThisFinalState(photonfs);
addProjection(vfs, "JetFS");
// Jets
FastJets jetpro(vfs, FastJets::ANTIKT, 0.4);
jetpro.useInvisibles();
addProjection(jetpro, "Jets");
_h_phbarrel_jetcentral_SS = bookHisto1D(1, 1, 1);
_h_phbarrel_jetmedium_SS = bookHisto1D(2, 1, 1);
_h_phbarrel_jetforward_SS = bookHisto1D(3, 1, 1);
_h_phbarrel_jetcentral_OS = bookHisto1D(4, 1, 1);
_h_phbarrel_jetmedium_OS = bookHisto1D(5, 1, 1);
_h_phbarrel_jetforward_OS = bookHisto1D(6, 1, 1);
}
int getEtaBin(double eta_w, int what) const {
double eta = fabs(eta_w);
int v_iter = 0;
if (what==0) {
for (v_iter=0; v_iter < (int)_eta_bins_ph.size()-1; v_iter++){
if (eta >= _eta_bins_ph.at(v_iter) && eta < _eta_bins_ph.at(v_iter+1))
break;
}
}
else if (what==1) {
for (v_iter=0; v_iter < (int)_eta_bins_jet.size()-1; v_iter++){
if (eta >= _eta_bins_jet.at(v_iter) && eta < _eta_bins_jet.at(v_iter+1))
break;
}
}
else {
for (v_iter=0; v_iter < (int)_eta_bins_areaoffset.size()-1; v_iter++){
if (eta >= _eta_bins_areaoffset.at(v_iter) && eta < _eta_bins_areaoffset.at(v_iter+1))
break;
}
}
return v_iter;
}
// Perform the per-event analysis
void analyze(const Event& event) {
const double weight = event.weight();
// Get the photon
const FinalState& photonfs = applyProjection<FinalState>(event, "LeadingPhoton");
if (photonfs.particles().size() < 1) vetoEvent;
const FourMomentum photon = photonfs.particles().front().momentum();
double eta_P = photon.eta();
double phi_P = photon.phi();
// Get the jet
Jets jets = applyProjection<FastJets>(event, "Jets").jetsByPt(20.0*GeV);
if (jets.size()==0) {
vetoEvent;
}
FourMomentum leadingJet = jets[0].momentum();
// Require jet separated from photon
if (deltaR(eta_P, phi_P, leadingJet.eta(), leadingJet.phi())<1.0) {
vetoEvent;
}
// Veto if leading jet is outside plotted rapidity regions
const double abs_y1 = fabs(leadingJet.rapidity());
if (abs_y1 > 4.4) {
vetoEvent;
}
// compute the median event energy density
const unsigned int skipnhardjets = 0;
_ptDensity.clear();
_sigma.clear();
_Njets.clear();
std::vector< std::vector<double> > ptDensities;
std::vector<double> emptyVec;
ptDensities.assign(_eta_bins_areaoffset.size()-1,emptyVec);
FastJets fast_jets = applyProjection<FastJets>(event, "KtJetsD05");
const fastjet::ClusterSequenceArea* clust_seq_area = fast_jets.clusterSeqArea();
foreach (const fastjet::PseudoJet& jet, fast_jets.pseudoJets(0.0*GeV)) {
double eta = fabs(jet.eta());
double pt = fabs(jet.perp());
double area = clust_seq_area->area(jet);
if (area > 10e-4 && fabs(eta)<_eta_bins_areaoffset[_eta_bins_areaoffset.size()-1]) {
ptDensities.at(getEtaBin(fabs(eta),2)).push_back(pt/area);
}
}
for (int b = 0; b<(int)_eta_bins_areaoffset.size()-1; b++) {
double median = 0.0;
double sigma = 0.0;
int Njets = 0;
if (ptDensities[b].size() > skipnhardjets) {
std::sort(ptDensities[b].begin(), ptDensities[b].end());
int nDens = ptDensities[b].size() - skipnhardjets;
if ( nDens%2 == 0 ) {
median = (ptDensities[b][nDens/2]+ptDensities[b][(nDens-2)/2])/2;
} else {
median = ptDensities[b][(nDens-1)/2];
}
sigma = ptDensities[b][(int)(.15865*nDens)];
Njets = nDens;
}
_ptDensity.push_back(median);
_sigma.push_back(sigma);
_Njets.push_back(Njets);
}
// compute photon isolation
// std EtCone
Particles fs = applyProjection<FinalState>(event, "JetFS").particles();
FourMomentum mom_in_EtCone;
float iso_dR = 0.4;
float cluster_eta_width = 0.25*7.0;
- float cluster_phi_width = (PI/128.)*5.0;
+ float cluster_phi_width = (PI/128.)*7.0;
foreach (const Particle& p, fs) {
// check if it's in the cone of .4
if (deltaR(eta_P, phi_P, p.eta(), p.phi()) >= iso_dR) continue;
// check if it's in the 5x7 central core
if (fabs(eta_P-p.eta()) < cluster_eta_width*0.5 &&
fabs(phi_P-p.phi()) < cluster_phi_width*0.5) continue;
mom_in_EtCone += p.momentum();
}
// now figure out the correction (area*density)
float EtCone_area = PI*iso_dR*iso_dR - cluster_eta_width*cluster_phi_width;
float correction = _ptDensity[getEtaBin(eta_P,2)]*EtCone_area;
// require photon to be isolated
if (mom_in_EtCone.Et()-correction > 4.0*GeV) vetoEvent;
int photon_jet_sign = sign( leadingJet.rapidity() * photon.rapidity() );
// Fill histos
float abs_jet_rapidity = fabs(leadingJet.rapidity());
float photon_pt = photon.pT()/GeV;
float abs_photon_eta = fabs(photon.eta());
if (abs_photon_eta<1.37) {
if (abs_jet_rapidity < 1.2) {
if (photon_jet_sign >= 1) {
_h_phbarrel_jetcentral_SS->fill(photon_pt, weight);
} else {
_h_phbarrel_jetcentral_OS->fill(photon_pt, weight);
}
} else if (abs_jet_rapidity < 2.8) {
if (photon_jet_sign >= 1) {
_h_phbarrel_jetmedium_SS->fill(photon_pt, weight);
} else {
_h_phbarrel_jetmedium_OS->fill(photon_pt, weight);
}
} else if (abs_jet_rapidity < 4.4) {
if (photon_jet_sign >= 1) {
_h_phbarrel_jetforward_SS->fill(photon_pt, weight);
} else {
_h_phbarrel_jetforward_OS->fill(photon_pt, weight);
}
}
}
}
/// Normalise histograms etc., after the run
void finalize() {
scale(_h_phbarrel_jetcentral_SS, crossSection()/sumOfWeights());
scale(_h_phbarrel_jetcentral_OS, crossSection()/sumOfWeights());
scale(_h_phbarrel_jetmedium_SS, crossSection()/sumOfWeights());
scale(_h_phbarrel_jetmedium_OS, crossSection()/sumOfWeights());
scale(_h_phbarrel_jetforward_SS, crossSection()/sumOfWeights());
scale(_h_phbarrel_jetforward_OS, crossSection()/sumOfWeights());
}
private:
Histo1DPtr _h_phbarrel_jetcentral_SS;
Histo1DPtr _h_phbarrel_jetmedium_SS;
Histo1DPtr _h_phbarrel_jetforward_SS;
Histo1DPtr _h_phbarrel_jetcentral_OS;
Histo1DPtr _h_phbarrel_jetmedium_OS;
Histo1DPtr _h_phbarrel_jetforward_OS;
fastjet::AreaDefinition* _area_def;
std::vector<float> _eta_bins_ph;
std::vector<float> _eta_bins_jet;
std::vector<float> _eta_bins_areaoffset;
std::vector<float> _ptDensity;
std::vector<float> _sigma;
std::vector<float> _Njets;
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(ATLAS_2012_I1093738);
}
diff --git a/src/Analyses/ATLAS_2012_I1199269.cc b/src/Analyses/ATLAS_2012_I1199269.cc
--- a/src/Analyses/ATLAS_2012_I1199269.cc
+++ b/src/Analyses/ATLAS_2012_I1199269.cc
@@ -1,193 +1,193 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
/// @todo These should be unnecessary from 2.2.0 onward when Rivet::Jet is rewritten
#include "fastjet/internal/base.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/AreaDefinition.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include "fastjet/PseudoJet.hh"
namespace Rivet {
/// @brief Measurement of isolated diphoton + X differential cross-sections
///
/// Inclusive isolated gamma gamma cross-sections, differential in M(gg), pT(gg),
/// dphi(gg), cos(theta*)_CS
///
/// @author Giovanni Marchiori
///
class ATLAS_2012_I1199269 : public Analysis {
public:
/// Constructor
ATLAS_2012_I1199269()
: Analysis("ATLAS_2012_I1199269")
{
_eta_bins_areaoffset += 0.0, 1.5, 3.0;
}
/// Book histograms and initialise projections before the run
void init() {
FinalState fs;
addProjection(fs, "FS");
/// @todo Clean-up when jets are better
FastJets fj(fs, FastJets::KT, 0.5);
_area_def = new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec());
fj.useJetArea(_area_def);
addProjection(fj, "KtJetsD05");
IdentifiedFinalState photonfs(Cuts::abseta < 2.37 && Cuts::pT > 22*GeV);
photonfs.acceptId(PID::PHOTON);
addProjection(photonfs, "Photon");
_h_M = bookHisto1D(1, 1, 1);
_h_pT = bookHisto1D(2, 1, 1);
_h_dPhi = bookHisto1D(3, 1, 1);
_h_cosThetaStar = bookHisto1D(4, 1, 1);
}
// int getBinIndex(double x, const vector<double>& binedges) const {
// /// @todo Add linear and log bin guessing, use binary split if sufficiently large, etc.
// for (size_t i = 0; i < _binedges.size()-1; ++i)
// if (inRange(x, binedges[i], binedges[i+1])) return (int) i;
// return -1;
// }
/// Perform the per-event analysis
void analyze(const Event& event) {
const double weight = event.weight();
/// Require at least 2 photons in final state
const Particles photons = applyProjection<IdentifiedFinalState>(event, "Photon").particlesByPt();
if (photons.size() < 2) vetoEvent;
/// Compute the median energy density
_ptDensity.clear();
_sigma.clear();
_Njets.clear();
/// @todo Nicer way to do this... in C++11?
vector<vector<double> > ptDensities;
vector<double> emptyVec;
ptDensities.assign(_eta_bins_areaoffset.size()-1, emptyVec);
// Get jets, and corresponding jet areas
const fastjet::ClusterSequenceArea* clust_seq_area = applyProjection<FastJets>(event, "KtJetsD05").clusterSeqArea();
foreach (const fastjet::PseudoJet& jet, applyProjection<FastJets>(event, "KtJetsD05").pseudoJets(0.0*GeV)) {
const double aeta = fabs(jet.eta());
const double pt = jet.perp();
const double area = clust_seq_area->area(jet);
if (area < 1e-3) continue;
const int ieta = binIndex(aeta, _eta_bins_areaoffset);
if (ieta != -1) ptDensities[ieta].push_back(pt/area);
}
// Compute median jet properties over the jets in the event
for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; ++b) {
double median = 0.0;
double sigma = 0.0;
int Njets = 0;
if (ptDensities[b].size() > 0) {
std::sort(ptDensities[b].begin(), ptDensities[b].end());
int nDens = ptDensities[b].size();
median = (nDens % 2 == 0) ? (ptDensities[b][nDens/2]+ptDensities[b][(nDens-2)/2])/2 : ptDensities[b][(nDens-1)/2];
sigma = ptDensities[b][(int)(.15865*nDens)];
Njets = nDens;
}
_ptDensity.push_back(median);
_sigma.push_back(sigma);
_Njets.push_back(Njets);
}
// Loop over photons and fill vector of isolated ones
Particles isolated_photons;
foreach (const Particle& photon, photons) {
/// Remove photons in ECAL crack region
if (inRange(photon.abseta(), 1.37, 1.52)) continue;
const double eta_P = photon.eta();
const double phi_P = photon.phi();
// Compute isolation via particles within an R=0.4 cone of the photon
Particles fs = applyProjection<FinalState>(event, "FS").particles();
FourMomentum mom_in_EtCone;
foreach (const Particle& p, fs) {
// Reject if not in cone
if (deltaR(photon.momentum(), p.momentum()) > 0.4) continue;
// Reject if in the 5x7 cell central core
if (fabs(eta_P - p.eta()) < 0.025 * 7 * 0.5 &&
- fabs(phi_P - p.phi()) < PI/128. * 5 * 0.5) continue;
+ fabs(phi_P - p.phi()) < PI/128. * 7 * 0.5) continue;
// Sum momentum
mom_in_EtCone += p.momentum();
}
// Now figure out the correction (area*density)
const double EtCone_area = PI*sqr(0.4) - (7*.025)*(5*PI/128.); // cone area - central core rectangle
const double correction = _ptDensity[binIndex(fabs(eta_P), _eta_bins_areaoffset)] * EtCone_area;
// Discard the photon if there is more than 4 GeV of cone activity
// NOTE: Shouldn't need to subtract photon itself (it's in the central core)
// NOTE: using expected cut at hadron/particle level, not at reco level
if (mom_in_EtCone.Et() - correction > 4*GeV) continue;
// Add isolated photon to list
isolated_photons.push_back(photon);
}
// Require at least two isolated photons and select leading pT pair
if (isolated_photons.size() < 2) vetoEvent;
sortByPt(isolated_photons);
FourMomentum y1 = isolated_photons[0].momentum();
FourMomentum y2 = isolated_photons[1].momentum();
// Leading photon should have pT > 25 GeV
if (y1.pT() < 25*GeV) vetoEvent;
// Require the two photons to be separated by dR > 0.4
if (deltaR(y1, y2) < 0.4) vetoEvent;
FourMomentum yy = y1 + y2;
const double Myy = yy.mass();
_h_M->fill(Myy/GeV, weight);
const double pTyy = yy.pT();
_h_pT->fill(pTyy/GeV, weight);
const double dPhiyy = mapAngle0ToPi(y1.phi() - y2.phi());
_h_dPhi->fill(dPhiyy, weight);
const double costhetayy = 2 * y1.pT() * y2.pT() * sinh(y1.eta() - y2.eta()) / Myy / add_quad(Myy, pTyy);
_h_cosThetaStar->fill(costhetayy, weight);
}
/// Normalise histograms etc., after the run
void finalize() {
scale(_h_M, crossSection()/sumOfWeights());
scale(_h_pT, crossSection()/sumOfWeights());
scale(_h_dPhi, crossSection()/sumOfWeights());
scale(_h_cosThetaStar, crossSection()/sumOfWeights());
}
private:
Histo1DPtr _h_M, _h_pT, _h_dPhi, _h_cosThetaStar;
fastjet::AreaDefinition* _area_def;
vector<double> _eta_bins;
vector<double> _eta_bins_areaoffset;
vector<double> _ptDensity;
vector<double> _sigma;
vector<double> _Njets;
};
DECLARE_RIVET_PLUGIN(ATLAS_2012_I1199269);
}
diff --git a/src/Analyses/ATLAS_2013_I1263495.cc b/src/Analyses/ATLAS_2013_I1263495.cc
--- a/src/Analyses/ATLAS_2013_I1263495.cc
+++ b/src/Analyses/ATLAS_2013_I1263495.cc
@@ -1,149 +1,149 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/LeadingParticlesFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
namespace Rivet {
/// Inclusive isolated prompt photon analysis with 2011 LHC data
class ATLAS_2013_I1263495 : public Analysis {
public:
/// Constructor
ATLAS_2013_I1263495()
: Analysis("ATLAS_2013_I1263495")
{
_eta_bins.push_back( 0.00);
_eta_bins.push_back( 1.37);
_eta_bins.push_back( 1.52);
_eta_bins.push_back( 2.37);
_eta_bins_areaoffset.push_back(0.0);
_eta_bins_areaoffset.push_back(1.5);
_eta_bins_areaoffset.push_back(3.0);
}
/// Book histograms and initialise projections before the run
void init() {
FinalState fs;
addProjection(fs, "FS");
// Consider the final state jets for the energy density calculation
FastJets fj(fs, FastJets::KT, 0.5);
/// @todo Memory leak! (Once per run, not serious)
_area_def = new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec()); ///< @todo FastJets should deep copy the pointer if possible
fj.useJetArea(_area_def);
addProjection(fj, "KtJetsD05");
// Consider the leading pt photon with |eta| < 2.37 and pT > 100 GeV
LeadingParticlesFinalState photonfs(FinalState(-2.37, 2.37, 100.0*GeV));
photonfs.addParticleId(PID::PHOTON);
addProjection(photonfs, "LeadingPhoton");
// Book the dsigma/dEt (in eta bins) histograms
for (size_t i = 0; i < _eta_bins.size() - 1; i++) {
if (fuzzyEquals(_eta_bins[i], 1.37)) continue; // skip this bin
_h_Et_photon[i] = bookHisto1D(1, 1, i+1);
}
// Book the dsigma/d|eta| histogram
_h_eta_photon = bookHisto1D(1,2,1);
}
/// Return eta bin for either dsigma/dET histogram (area_eta=false) or energy density correction (area_eta=true)
size_t _getEtaBin(double eta_w, bool area_eta) const {
const double eta = fabs(eta_w);
if (!area_eta) {
return binIndex(eta, _eta_bins);
} else {
return binIndex(eta, _eta_bins_areaoffset);
}
}
/// Perform the per-event analysis
void analyze(const Event& event) {
// Retrieve leading photon
Particles photons = applyProjection<LeadingParticlesFinalState>(event, "LeadingPhoton").particles();
if (photons.size() != 1) vetoEvent;
const Particle& leadingPhoton = photons[0];
// Veto events with photon in ECAL crack
if (inRange(leadingPhoton.abseta(), 1.37, 1.52)) vetoEvent;
// Compute isolation energy in cone of radius .4 around photon (all particles)
FourMomentum mom_in_EtCone;
Particles fs = applyProjection<FinalState>(event, "FS").particles();
foreach (const Particle& p, fs) {
// Check if it's outside the cone of 0.4
if (deltaR(leadingPhoton, p) >= 0.4) continue;
// Don't count particles in the 5x7 central core
if (deltaEta(leadingPhoton, p) < .025*7.0*0.5 &&
- deltaPhi(leadingPhoton, p) < (PI/128.)*5.0*0.5) continue;
+ deltaPhi(leadingPhoton, p) < (PI/128.)*7.0*0.5) continue;
// Increment isolation energy
mom_in_EtCone += p.momentum();
}
// Get the area-filtered jet inputs for computing median energy density, etc.
vector<double> ptDensity, ptSigma, nJets;
vector< vector<double> > ptDensities(_eta_bins_areaoffset.size()-1);
FastJets fast_jets =applyProjection<FastJets>(event, "KtJetsD05");
const fastjet::ClusterSequenceArea* clust_seq_area = fast_jets.clusterSeqArea();
foreach (const Jet& jet, fast_jets.jets()) {
const double area = clust_seq_area->area(jet);
if (area > 10e-4 && jet.abseta() < _eta_bins_areaoffset.back())
ptDensities.at( _getEtaBin(jet.abseta(), true) ).push_back(jet.pT()/area);
}
// Compute the median energy density, etc.
for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; b++) {
const int njets = ptDensities[b].size();
const double ptmedian = (njets > 0) ? median(ptDensities[b]) : 0;
const double ptsigma = (njets > 0) ? ptDensities[b][(size_t)(0.15865*njets)] : 0;
nJets.push_back(njets);
ptDensity.push_back(ptmedian);
ptSigma.push_back(ptsigma);
}
// Compute the isolation energy correction (cone area*energy density)
const double etCone_area = PI*sqr(0.4) - (7.0*.025)*(5.0*PI/128.);
const double correction = ptDensity[_getEtaBin(leadingPhoton.abseta(), true)]*etCone_area;
// Apply isolation cut on area-corrected value
if (mom_in_EtCone.Et() - correction > 7*GeV) vetoEvent;
// Fill histograms
const size_t eta_bin = _getEtaBin(leadingPhoton.abseta(), false);
_h_Et_photon[eta_bin]->fill(leadingPhoton.Et(), event.weight());
_h_eta_photon->fill(leadingPhoton.abseta(), event.weight());
}
/// Normalise histograms etc., after the run
void finalize() {
for (size_t i = 0; i < _eta_bins.size()-1; i++) {
if (fuzzyEquals(_eta_bins[i], 1.37)) continue;
scale(_h_Et_photon[i], crossSection()/picobarn/sumOfWeights());
}
scale(_h_eta_photon, crossSection()/picobarn/sumOfWeights());
}
private:
Histo1DPtr _h_Et_photon[3];
Histo1DPtr _h_eta_photon;
fastjet::AreaDefinition* _area_def;
vector<double> _eta_bins, _eta_bins_areaoffset;
};
DECLARE_RIVET_PLUGIN(ATLAS_2013_I1263495);
}
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Tue, Nov 19, 8:55 PM (1 d, 2 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806180
Default Alt Text
(41 KB)
Attached To
rRIVETHG rivethg
Event Timeline
Log In to Comment