Page MenuHomeHEPForge
Files F7879209

No OneTemporary
Actions

View Options

diff --git a/analyses/pluginALICE/ALICE_2010_S8624100.cc b/analyses/pluginALICE/ALICE_2010_S8624100.cc
--- a/analyses/pluginALICE/ALICE_2010_S8624100.cc
+++ b/analyses/pluginALICE/ALICE_2010_S8624100.cc
@@ -1,94 +1,92 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
namespace Rivet {
class ALICE_2010_S8624100 : public Analysis {
public:
/// @name Constructors etc.
//@{
/// Constructor
ALICE_2010_S8624100()
: Analysis("ALICE_2010_S8624100")
{ }
//@}
public:
/// @name Analysis methods
//@{
/// Book histograms and initialise projections before the run
void init() {
ChargedFinalState cfs05(-0.5, 0.5);
ChargedFinalState cfs10(-1.0, 1.0);
ChargedFinalState cfs13(-1.3, 1.3);
declare(cfs05, "CFS05");
declare(cfs10, "CFS10");
declare(cfs13, "CFS13");
if (fuzzyEquals(sqrtS()/GeV, 900, 1E-3)) {
book(_h_dN_dNch_05 ,11, 1, 1);
book(_h_dN_dNch_10 ,12, 1, 1);
book(_h_dN_dNch_13 ,13, 1, 1);
} else if (fuzzyEquals(sqrtS()/GeV, 2360, 1E-3)) {
book(_h_dN_dNch_05 ,17, 1, 1);
book(_h_dN_dNch_10 ,18, 1, 1);
book(_h_dN_dNch_13 ,19, 1, 1);
}
}
/// Perform the per-event analysis
void analyze(const Event& event) {
- const double weight = 1.0;
-
- const ChargedFinalState& charged_05 = apply<ChargedFinalState>(event, "CFS05");
+ const ChargedFinalState& charged_05 = apply<ChargedFinalState>(event, "CFS05");
const ChargedFinalState& charged_10 = apply<ChargedFinalState>(event, "CFS10");
const ChargedFinalState& charged_13 = apply<ChargedFinalState>(event, "CFS13");
- _h_dN_dNch_05->fill(charged_05.size(), weight);
- _h_dN_dNch_10->fill(charged_10.size(), weight);
- _h_dN_dNch_13->fill(charged_13.size(), weight);
+ _h_dN_dNch_05->fill(charged_05.size());
+ _h_dN_dNch_10->fill(charged_10.size());
+ _h_dN_dNch_13->fill(charged_13.size());
}
/// Normalise histograms etc., after the run
void finalize() {
normalize(_h_dN_dNch_05);
normalize(_h_dN_dNch_10);
normalize(_h_dN_dNch_13);
}
//@}
private:
/// @name Histograms
//@{
Histo1DPtr _h_dN_dNch_05;
Histo1DPtr _h_dN_dNch_10;
Histo1DPtr _h_dN_dNch_13;
//@}
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(ALICE_2010_S8624100);
}
diff --git a/analyses/pluginALICE/ALICE_2010_S8625980.cc b/analyses/pluginALICE/ALICE_2010_S8625980.cc
--- a/analyses/pluginALICE/ALICE_2010_S8625980.cc
+++ b/analyses/pluginALICE/ALICE_2010_S8625980.cc
@@ -1,99 +1,97 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
namespace Rivet {
class ALICE_2010_S8625980 : public Analysis {
public:
/// @name Constructors etc.
//@{
/// Constructor
ALICE_2010_S8625980()
- : Analysis("ALICE_2010_S8625980"),
- _Nevt_after_cuts(0.0)
+ : Analysis("ALICE_2010_S8625980")
{ }
//@}
public:
/// @name Analysis methods
//@{
/// Book histograms and initialise projections before the run
void init() {
ChargedFinalState cfs(-1.0, 1.0);
declare(cfs, "CFS");
if (fuzzyEquals(sqrtS()/GeV, 900, 1E-3)) {
book(_h_dN_deta ,4, 1, 1);
} else if (fuzzyEquals(sqrtS()/GeV, 2360, 1E-3)) {
book(_h_dN_deta ,5, 1, 1);
} else if (fuzzyEquals(sqrtS()/GeV, 7000, 1E-3)) {
book(_h_dN_deta ,6, 1, 1);
book(_h_dN_dNch ,3, 1, 1);
}
+ book(_Nevt_after_cuts, "Nevt_after_cuts");
}
/// Perform the per-event analysis
void analyze(const Event& event) {
- const double weight = 1.0;
-
const ChargedFinalState& charged = apply<ChargedFinalState>(event, "CFS");
if (charged.size() < 1) {
vetoEvent;
}
- _Nevt_after_cuts += weight;
+ _Nevt_after_cuts->fill();
foreach (const Particle& p, charged.particles()) {
const double eta = p.eta();
- _h_dN_deta->fill(eta, weight);
+ _h_dN_deta->fill(eta);
}
if (fuzzyEquals(sqrtS()/GeV, 7000, 1E-3)) {
- _h_dN_dNch->fill(charged.size(), weight);
+ _h_dN_dNch->fill(charged.size());
}
}
/// Normalise histograms etc., after the run
void finalize() {
if (fuzzyEquals(sqrtS()/GeV, 7000, 1E-3)) {
normalize(_h_dN_dNch);
}
scale(_h_dN_deta, 1.0/_Nevt_after_cuts);
}
//@}
private:
/// @name Histograms
//@{
Histo1DPtr _h_dN_deta;
Histo1DPtr _h_dN_dNch;
- double _Nevt_after_cuts;
+ CounterPtr _Nevt_after_cuts;
//@}
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(ALICE_2010_S8625980);
}
diff --git a/analyses/pluginALICE/ALICE_2010_S8706239.cc b/analyses/pluginALICE/ALICE_2010_S8706239.cc
--- a/analyses/pluginALICE/ALICE_2010_S8706239.cc
+++ b/analyses/pluginALICE/ALICE_2010_S8706239.cc
@@ -1,101 +1,100 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
namespace Rivet {
class ALICE_2010_S8706239 : public Analysis {
public:
/// @name Constructors etc.
//@{
/// Constructor
ALICE_2010_S8706239()
- : Analysis("ALICE_2010_S8706239"),
- _Nevt_after_cuts(0.0)
+ : Analysis("ALICE_2010_S8706239")
{ }
//@}
public:
/// @name Analysis methods
//@{
/// Book histograms and initialise projections before the run
void init() {
ChargedFinalState cfs(-0.8, 0.8, 0.15);
declare(cfs, "CFS");
book(_h_pT ,4, 1, 1);
book(_h_pT_Nch_015 ,11, 1, 1);
book(_h_pT_Nch_05 ,12, 1, 1);
+ book(_Nevt_after_cuts,"Nevt_after_cuts");
+
}
/// Perform the per-event analysis
void analyze(const Event& event) {
- const double weight = 1.0;
-
const ChargedFinalState& charged = apply<ChargedFinalState>(event, "CFS");
- _Nevt_after_cuts += weight;
+ _Nevt_after_cuts->fill();
// Get number of particles that fulfill certain pT requirements
int Nch_015 = 0;
int Nch_05 = 0;
foreach (const Particle& p, charged.particles()) {
double pT = p.pT()/GeV;
if (pT < 4.0) Nch_015++;
if (pT > 0.5 && pT < 4.0) Nch_05++;
}
// Now we can fill histograms
foreach (const Particle& p, charged.particles()) {
double pT = p.pT()/GeV;
- if (pT < 4.0) _h_pT_Nch_015 ->fill(Nch_015, pT, weight);
- if (pT > 0.5 && pT < 4.0) _h_pT_Nch_05 ->fill(Nch_05, pT, weight);
+ if (pT < 4.0) _h_pT_Nch_015 ->fill(Nch_015, pT);
+ if (pT > 0.5 && pT < 4.0) _h_pT_Nch_05 ->fill(Nch_05, pT);
// To get the Yield, fill appropriate weight 1/(2PI * pT * d eta)
- _h_pT->fill(pT, weight /(TWOPI*pT*1.6) );
+ _h_pT->fill(pT, 1.0 /(TWOPI*pT*1.6) );
}
}
/// Normalise histograms etc., after the run
void finalize() {
scale(_h_pT, 1.0/_Nevt_after_cuts);
}
//@}
private:
/// @name Histograms
//@{
Histo1DPtr _h_pT;
Profile1DPtr _h_pT_Nch_015 ;
Profile1DPtr _h_pT_Nch_05 ;
- double _Nevt_after_cuts;
+ CounterPtr _Nevt_after_cuts;
//@}
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(ALICE_2010_S8706239);
}
diff --git a/analyses/pluginALICE/ALICE_2011_S8909580.cc b/analyses/pluginALICE/ALICE_2011_S8909580.cc
--- a/analyses/pluginALICE/ALICE_2011_S8909580.cc
+++ b/analyses/pluginALICE/ALICE_2011_S8909580.cc
@@ -1,103 +1,102 @@
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/UnstableFinalState.hh"
namespace Rivet {
class ALICE_2011_S8909580 : public Analysis {
public:
ALICE_2011_S8909580()
: Analysis("ALICE_2011_S8909580")
{}
public:
void init() {
const UnstableFinalState ufs(Cuts::abseta < 15);
declare(ufs, "UFS");
book(_histPtK0s ,1, 1, 1);
book(_histPtLambda ,2, 1, 1);
book(_histPtAntiLambda ,3, 1, 1);
book(_histPtXi ,4, 1, 1);
book(_histPtPhi ,5, 1, 1);
book(_temp_h_Lambdas ,"TMP/h_Lambdas", refData(6, 1, 1));
book(_temp_h_Kzeros ,"TMP/h_Kzeros", refData(6, 1, 1));
book(_h_LamKzero ,6, 1, 1);
}
void analyze(const Event& event) {
- const double weight = 1.0;
const UnstableFinalState& ufs = apply<UnstableFinalState>(event, "UFS");
foreach (const Particle& p, ufs.particles()) {
const double absrap = p.absrap();
const double pT = p.pT()/GeV;
if (absrap < 0.8) {
switch(p.pid()) {
case 3312:
case -3312:
if ( !( p.hasAncestor(3334) || p.hasAncestor(-3334) ) ) {
- _histPtXi->fill(pT, weight);
+ _histPtXi->fill(pT);
}
break;
if (absrap < 0.75) {
case 310:
- _histPtK0s->fill(pT, weight);
- _temp_h_Kzeros->fill(pT, 2*weight);
+ _histPtK0s->fill(pT);
+ _temp_h_Kzeros->fill(pT, 2);
break;
case 3122:
if ( !( p.hasAncestor(3322) || p.hasAncestor(-3322) ||
p.hasAncestor(3312) || p.hasAncestor(-3312) ||
p.hasAncestor(3334) || p.hasAncestor(-3334) ) ) {
- _histPtLambda->fill(pT, weight);
- _temp_h_Lambdas->fill(pT, weight);
+ _histPtLambda->fill(pT);
+ _temp_h_Lambdas->fill(pT);
}
break;
case -3122:
if ( !( p.hasAncestor(3322) || p.hasAncestor(-3322) ||
p.hasAncestor(3312) || p.hasAncestor(-3312) ||
p.hasAncestor(3334) || p.hasAncestor(-3334) ) ) {
- _histPtAntiLambda->fill(pT, weight);
- _temp_h_Lambdas->fill(pT, weight);
+ _histPtAntiLambda->fill(pT);
+ _temp_h_Lambdas->fill(pT);
}
break;
}
if (absrap<0.6) {
case 333:
- _histPtPhi->fill(pT, weight);
+ _histPtPhi->fill(pT);
break;
}
}
}
}
}
void finalize() {
scale(_histPtK0s, 1./(1.5*sumOfWeights()));
scale(_histPtLambda, 1./(1.5*sumOfWeights()));
scale(_histPtAntiLambda, 1./(1.5*sumOfWeights()));
scale(_histPtXi, 1./(1.6*sumOfWeights()));
scale(_histPtPhi, 1./(1.2*sumOfWeights()));
divide(_temp_h_Lambdas, _temp_h_Kzeros, _h_LamKzero);
}
private:
Histo1DPtr _histPtK0s, _histPtLambda, _histPtAntiLambda, _histPtXi, _histPtPhi;
Histo1DPtr _temp_h_Lambdas, _temp_h_Kzeros;
Scatter2DPtr _h_LamKzero;
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(ALICE_2011_S8909580);
}
diff --git a/analyses/pluginALICE/ALICE_2011_S8945144.cc b/analyses/pluginALICE/ALICE_2011_S8945144.cc
--- a/analyses/pluginALICE/ALICE_2011_S8945144.cc
+++ b/analyses/pluginALICE/ALICE_2011_S8945144.cc
@@ -1,105 +1,104 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
namespace Rivet {
class ALICE_2011_S8945144 : public Analysis {
public:
ALICE_2011_S8945144()
: Analysis("ALICE_2011_S8945144")
{}
public:
void init() {
const ChargedFinalState cfs(-15, 15);
declare(cfs, "CFS");
book(_histPtPions ,"d01-x01-y01");
book(_histPtAntiPions ,"d01-x01-y02");
book(_histPtKaons ,"d02-x01-y01");
book(_histPtAntiKaons ,"d02-x01-y02");
book(_histPtProtons ,"d03-x01-y01");
book(_histPtAntiProtons ,"d03-x01-y02");
book(_histAveragePt ,"d04-x01-y01");
}
void analyze(const Event& event) {
- const double weight = 1.0;
const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
foreach (const Particle& p, cfs.particles()) {
if(p.absrap()<0.5) {
switch (p.pid()) {
case 211:
- _histPtPions->fill(p.pT()/GeV, weight);
- _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV, weight);
+ _histPtPions->fill(p.pT()/GeV);
+ _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV);
break;
case -211:
- _histPtAntiPions->fill(p.pT()/GeV, weight);
- _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV, weight);
+ _histPtAntiPions->fill(p.pT()/GeV);
+ _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV);
break;
case 2212:
if ( !(p.hasAncestor(3322) || // Xi0
p.hasAncestor(3122) || p.hasAncestor(-3122) || // Lambda
p.hasAncestor(3222) || p.hasAncestor(-3222) || // Sigma+/-
p.hasAncestor(3312) || p.hasAncestor(-3312) ) ) { // Xi-/+
- _histPtProtons->fill(p.pT()/GeV, weight);
- _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV, weight);
+ _histPtProtons->fill(p.pT()/GeV);
+ _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV);
}
break;
case -2212:
if ( !(p.hasAncestor(3322) || // Xi0
p.hasAncestor(3122) || p.hasAncestor(-3122) || // Lambda
p.hasAncestor(3222) || p.hasAncestor(-3222) || // Sigma+/-
p.hasAncestor(3312) || p.hasAncestor(-3312) ) ) { // Xi-/+
- _histPtAntiProtons->fill(p.pT()/GeV, weight);
- _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV, weight);
+ _histPtAntiProtons->fill(p.pT()/GeV);
+ _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV);
}
break;
case 321:
- _histPtKaons->fill(p.pT()/GeV, weight);
- _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV, weight);
+ _histPtKaons->fill(p.pT()/GeV);
+ _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV);
break;
case -321:
- _histPtAntiKaons->fill(p.pT()/GeV, weight);
- _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV, weight);
+ _histPtAntiKaons->fill(p.pT()/GeV);
+ _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV);
break;
}
}
}
}
void finalize() {
scale(_histPtPions, 1./sumOfWeights());
scale(_histPtProtons, 1./sumOfWeights());
scale(_histPtKaons, 1./sumOfWeights());
scale(_histPtAntiPions, 1./sumOfWeights());
scale(_histPtAntiProtons, 1./sumOfWeights());
scale(_histPtAntiKaons, 1./sumOfWeights());
}
private:
Histo1DPtr _histPtPions;
Histo1DPtr _histPtProtons;
Histo1DPtr _histPtKaons;
Histo1DPtr _histPtAntiPions;
Histo1DPtr _histPtAntiProtons;
Histo1DPtr _histPtAntiKaons;
Profile1DPtr _histAveragePt;
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(ALICE_2011_S8945144);
}
diff --git a/analyses/pluginALICE/ALICE_2012_I1181770.cc b/analyses/pluginALICE/ALICE_2012_I1181770.cc
--- a/analyses/pluginALICE/ALICE_2012_I1181770.cc
+++ b/analyses/pluginALICE/ALICE_2012_I1181770.cc
@@ -1,116 +1,114 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
namespace Rivet {
class ALICE_2012_I1181770 : public Analysis {
public:
ALICE_2012_I1181770()
: Analysis("ALICE_2012_I1181770")
{ }
void init() {
// Projection setup
declare(ChargedFinalState(), "CFS");
// Book (energy-specific) histograms
int isqrts = -1;
if (fuzzyEquals(sqrtS()/GeV, 900, 1E-3)) isqrts = 1;
else if (fuzzyEquals(sqrtS()/GeV, 2760, 1E-3)) isqrts = 2;
else if (fuzzyEquals(sqrtS()/GeV, 7000, 1E-3)) isqrts = 3;
assert(isqrts > 0);
book(_h_frac_sd_inel, 1, 1, isqrts);
book(_h_frac_dd_inel, 2, 1, isqrts);
book(_h_xsec_sd , 3, 1, isqrts);
book(_h_xsec_dd , 4, 1, isqrts);
book(_h_xsec_inel , 5, 1, isqrts);
}
void analyze(const Event& event) {
const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
if (cfs.size() < 2) vetoEvent; // need at least two particles to calculate gaps
- const double weight = 1.0;
-
// Fill INEL plots for each event
- _h_xsec_inel->fill(sqrtS()/GeV, weight);
+ _h_xsec_inel->fill(sqrtS()/GeV);
// Identify particles with most positive/most negative rapidities
const Particles particlesByRap = cfs.particles(cmpMomByRap);
const Particle pslowest = particlesByRap.front();
const Particle pfastest = particlesByRap.back();
// Find gap sizes
const Particles particlesByEta = cfs.particles(cmpMomByEta); // sorted from minus to plus
const size_t num_particles = particlesByEta.size();
vector<double> gaps;
for (size_t ip = 1; ip < num_particles; ++ip) {
const Particle& p1 = particlesByEta[ip-1];
const Particle& p2 = particlesByEta[ip];
const double gap = p2.eta() - p1.eta();
assert(gap >= 0);
gaps.push_back(gap);
}
// First, last, and largest gaps
const double gapmax = *max_element(gaps.begin(), gaps.end());
const double gapbwd = gaps.front();
const double gapfwd = gaps.back();
// Mx calculation
FourMomentum p4lead;
if (pslowest.pid() == PID::PROTON && pfastest.pid() == PID::PROTON) {
p4lead = (fabs(pslowest.rapidity()) > fabs(pfastest.rapidity())) ? pslowest.momentum() : pfastest.momentum();
} else if (pslowest.pid() == PID::PROTON) {
p4lead = pslowest.momentum();
} else if (pfastest.pid() == PID::PROTON) {
p4lead = pfastest.momentum();
}
const double Mx = sqrt( (sqrtS()-p4lead.E()-p4lead.p3().mod()) * (sqrtS()-p4lead.E()+p4lead.p3().mod()) );
// Fill SD (and escape) if Mx is sufficiently low
if (Mx < 200*GeV) {
- _h_xsec_sd->fill(sqrtS()/GeV, weight);
+ _h_xsec_sd->fill(sqrtS()/GeV);
return;
}
// Also remove SD-like events in NSD events
if (fuzzyEquals(gapbwd, gapmax) || fuzzyEquals(gapfwd, gapmax)) vetoEvent;
// Fill DD plots
- if (gapmax > 3) _h_xsec_dd->fill(sqrtS()/GeV, weight);
+ if (gapmax > 3) _h_xsec_dd->fill(sqrtS()/GeV);
}
void finalize() {
// get the ratio plots: SD/inel, DD/inel
divide(_h_xsec_sd , _h_xsec_inel, _h_frac_sd_inel);
divide(_h_xsec_sd , _h_xsec_inel, _h_frac_dd_inel);
const double scaling = crossSection()/millibarn/sumOfWeights();
scale(_h_xsec_sd, scaling);
scale(_h_xsec_dd, scaling);
scale(_h_xsec_inel, scaling);
}
private:
Scatter2DPtr _h_frac_sd_inel;
Scatter2DPtr _h_frac_dd_inel;
Histo1DPtr _h_xsec_sd;
Histo1DPtr _h_xsec_dd;
Histo1DPtr _h_xsec_inel;
};
// Hook for the plugin system
DECLARE_RIVET_PLUGIN(ALICE_2012_I1181770);
}
diff --git a/analyses/pluginALICE/ALICE_2014_I1300380.cc b/analyses/pluginALICE/ALICE_2014_I1300380.cc
--- a/analyses/pluginALICE/ALICE_2014_I1300380.cc
+++ b/analyses/pluginALICE/ALICE_2014_I1300380.cc
@@ -1,120 +1,119 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/UnstableFinalState.hh"
namespace Rivet {
class ALICE_2014_I1300380 : public Analysis {
public:
ALICE_2014_I1300380()
: Analysis("ALICE_2014_I1300380")
{}
public:
void init() {
const UnstableFinalState cfs(Cuts::absrap<0.5);
declare(cfs, "CFS");
// Plots from the paper
book(_histPtSigmaStarPlus ,"d01-x01-y01"); // Sigma*+
book(_histPtSigmaStarMinus ,"d01-x01-y02"); // Sigma*-
book(_histPtSigmaStarPlusAnti ,"d01-x01-y03"); // anti Sigma*-
book(_histPtSigmaStarMinusAnti ,"d01-x01-y04"); // anti Sigma*+
book(_histPtXiStar ,"d02-x01-y01"); // 0.5 * (xi star + anti xi star)
book(_histAveragePt ,"d03-x01-y01"); // <pT> profile
}
void analyze(const Event& event) {
- const double weight = 1.0;
const UnstableFinalState& cfs = apply<UnstableFinalState>(event, "CFS");
foreach (const Particle& p, cfs.particles()) {
// protections against mc generators decaying long-lived particles
if ( !(p.hasAncestor(310) || p.hasAncestor(-310) || // K0s
p.hasAncestor(130) || p.hasAncestor(-130) || // K0l
p.hasAncestor(3322) || p.hasAncestor(-3322) || // Xi0
p.hasAncestor(3122) || p.hasAncestor(-3122) || // Lambda
p.hasAncestor(3222) || p.hasAncestor(-3222) || // Sigma+/-
p.hasAncestor(3312) || p.hasAncestor(-3312) || // Xi-/+
p.hasAncestor(3334) || p.hasAncestor(-3334) )) // Omega-/+
{
int aid = abs(p.pdgId());
if (aid == 211 || // pi+
aid == 321 || // K+
aid == 313 || // K*(892)0
aid == 2212 || // proton
aid == 333 ) { // phi(1020)
- _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV, weight);
+ _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV);
}
} // end if "rejection of long-lived particles"
switch (p.pdgId()) {
case 3224:
- _histPtSigmaStarPlus->fill(p.pT()/GeV, weight);
- _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV, weight);
+ _histPtSigmaStarPlus->fill(p.pT()/GeV);
+ _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV);
break;
case -3224:
- _histPtSigmaStarPlusAnti->fill(p.pT()/GeV, weight);
- _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV, weight);
+ _histPtSigmaStarPlusAnti->fill(p.pT()/GeV);
+ _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV);
break;
case 3114:
- _histPtSigmaStarMinus->fill(p.pT()/GeV, weight);
- _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV, weight);
+ _histPtSigmaStarMinus->fill(p.pT()/GeV);
+ _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV);
break;
case -3114:
- _histPtSigmaStarMinusAnti->fill(p.pT()/GeV, weight);
- _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV, weight);
+ _histPtSigmaStarMinusAnti->fill(p.pT()/GeV);
+ _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV);
break;
case 3324:
- _histPtXiStar->fill(p.pT()/GeV, weight);
- _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV, weight);
+ _histPtXiStar->fill(p.pT()/GeV);
+ _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV);
break;
case -3324:
- _histPtXiStar->fill(p.pT()/GeV, weight);
- _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV, weight);
+ _histPtXiStar->fill(p.pT()/GeV);
+ _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV);
break;
case 3312:
- _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV, weight);
+ _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV);
break;
case -3312:
- _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV, weight);
+ _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV);
break;
case 3334:
- _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV, weight);
+ _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV);
break;
case -3334:
- _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV, weight);
+ _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV);
break;
}
}
}
void finalize() {
scale(_histPtSigmaStarPlus, 1./sumOfWeights());
scale(_histPtSigmaStarPlusAnti, 1./sumOfWeights());
scale(_histPtSigmaStarMinus, 1./sumOfWeights());
scale(_histPtSigmaStarMinusAnti, 1./sumOfWeights());
scale(_histPtXiStar, 1./sumOfWeights()/ 2.);
}
private:
// plots from the paper
Histo1DPtr _histPtSigmaStarPlus;
Histo1DPtr _histPtSigmaStarPlusAnti;
Histo1DPtr _histPtSigmaStarMinus;
Histo1DPtr _histPtSigmaStarMinusAnti;
Histo1DPtr _histPtXiStar;
Profile1DPtr _histAveragePt;
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(ALICE_2014_I1300380);
}
diff --git a/analyses/pluginALICE/ALICE_2015_I1357424.cc b/analyses/pluginALICE/ALICE_2015_I1357424.cc
--- a/analyses/pluginALICE/ALICE_2015_I1357424.cc
+++ b/analyses/pluginALICE/ALICE_2015_I1357424.cc
@@ -1,99 +1,98 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
namespace Rivet {
class ALICE_2015_I1357424 : public Analysis {
public:
ALICE_2015_I1357424()
: Analysis("ALICE_2015_I1357424")
{}
public:
void init() {
const ChargedFinalState cfs(Cuts::absrap<0.5);
declare(cfs, "CFS");
//
// plots from the paper
book(_histPtPions ,"d01-x01-y01"); // pions
book(_histPtKaons ,"d01-x01-y02"); // kaons
book(_histPtProtons ,"d01-x01-y03"); // protons
book(_histPtKtoPi ,"d02-x01-y01"); // K to pi ratio
book(_histPtPtoPi ,"d03-x01-y01"); // p to pi ratio
//
// temp histos for ratios
book(_histPtPionsR1 ,"TMP/pT_pi1", refData(2, 1, 1)); // pi histo compatible with more restricted kaon binning
book(_histPtPionsR2 ,"TMP/pT_pi2", refData(3, 1, 1)); // pi histo compatible with more restricted proton binning
book(_histPtKaonsR ,"TMP/pT_K", refData(2, 1, 1)); // K histo with more restricted binning
book(_histPtProtonsR ,"TMP/pT_p", refData(3, 1, 1)); // p histo with more restricted binning
}
void analyze(const Event& event) {
- const double weight = 1.0;
const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
foreach (const Particle& p, cfs.particles()) {
// protections against mc generators decaying long-lived particles
if ( !(p.hasAncestor(310) || p.hasAncestor(-310) || // K0s
p.hasAncestor(130) || p.hasAncestor(-130) || // K0l
p.hasAncestor(3322) || p.hasAncestor(-3322) || // Xi0
p.hasAncestor(3122) || p.hasAncestor(-3122) || // Lambda
p.hasAncestor(3222) || p.hasAncestor(-3222) || // Sigma+/-
p.hasAncestor(3312) || p.hasAncestor(-3312) || // Xi-/+
p.hasAncestor(3334) || p.hasAncestor(-3334) )) // Omega-/+
{
switch (abs(p.pid())) {
case 211: // pi+
- _histPtPions->fill(p.pT()/GeV, weight);
- _histPtPionsR1->fill(p.pT()/GeV, weight);
- _histPtPionsR2->fill(p.pT()/GeV, weight);
+ _histPtPions->fill(p.pT()/GeV);
+ _histPtPionsR1->fill(p.pT()/GeV);
+ _histPtPionsR2->fill(p.pT()/GeV);
break;
case 2212: // proton
- _histPtProtons->fill(p.pT()/GeV, weight);
- _histPtProtonsR->fill(p.pT()/GeV, weight);
+ _histPtProtons->fill(p.pT()/GeV);
+ _histPtProtonsR->fill(p.pT()/GeV);
break;
case 321: // K+
- _histPtKaons->fill(p.pT()/GeV, weight);
- _histPtKaonsR->fill(p.pT()/GeV, weight);
+ _histPtKaons->fill(p.pT()/GeV);
+ _histPtKaonsR->fill(p.pT()/GeV);
break;
} // particle switch
} // primary pi, K, p only
} // particle loop
}
void finalize() {
divide(_histPtKaonsR, _histPtPionsR1, _histPtKtoPi);
divide(_histPtProtonsR, _histPtPionsR2, _histPtPtoPi);
scale(_histPtPions, 1./sumOfWeights());
scale(_histPtProtons, 1./sumOfWeights());
scale(_histPtKaons, 1./sumOfWeights());
}
private:
Histo1DPtr _histPtPions;
Histo1DPtr _histPtProtons;
Histo1DPtr _histPtKaons;
Histo1DPtr _histPtPionsR1;
Histo1DPtr _histPtPionsR2;
Histo1DPtr _histPtProtonsR;
Histo1DPtr _histPtKaonsR;
Scatter2DPtr _histPtKtoPi;
Scatter2DPtr _histPtPtoPi;
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(ALICE_2015_I1357424);
}
diff --git a/analyses/pluginATLAS/ATLAS_2012_I1124167.cc b/analyses/pluginATLAS/ATLAS_2012_I1124167..cc.needstransversesphericity
rename from analyses/pluginATLAS/ATLAS_2012_I1124167.cc
rename to analyses/pluginATLAS/ATLAS_2012_I1124167..cc.needstransversesphericity
diff --git a/analyses/pluginLHCb/LHCB_2011_I917009.cc b/analyses/pluginLHCb/LHCB_2011_I917009.cc
--- a/analyses/pluginLHCb/LHCB_2011_I917009.cc
+++ b/analyses/pluginLHCb/LHCB_2011_I917009.cc
@@ -1,323 +1,334 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/UnstableFinalState.hh"
namespace Rivet {
class LHCB_2011_I917009 : public Analysis {
public:
/// @name Constructors etc.
//@{
/// Constructor
LHCB_2011_I917009()
: Analysis("LHCB_2011_I917009"),
rap_beam(0.0), pt_min(0.0),
pt1_edge(0.65), pt2_edge(1.0),
pt3_edge(2.5), rap_min(2.),
rap_max(0.0), dsShift(0)
{ }
//@}
public:
/// @name Analysis methods
//@{
/// Book histograms and initialise projections before the run
void init() {
int y_nbins = 4;
fillMap(partLftMap);
if (fuzzyEquals(sqrtS(), 0.9*TeV)) {
rap_beam = 6.87;
rap_max = 4.;
pt_min = 0.25;
} else if (fuzzyEquals(sqrtS(), 7*TeV)) {
rap_beam = 8.92;
rap_max = 4.5;
pt_min = 0.15;
y_nbins = 5;
dsShift = 8;
} else {
MSG_ERROR("Incompatible beam energy!");
}
// Create the sets of temporary histograms that will be used to make the ratios in the finalize()
- for (size_t i = 0; i < 12; ++i) _tmphistos[i] = YODA::Histo1D(y_nbins, rap_min, rap_max);
- for (size_t i = 12; i < 15; ++i) _tmphistos[i] = YODA::Histo1D(refData(dsShift+5, 1, 1));
- for (size_t i = 15; i < 18; ++i) _tmphistos[i] = YODA::Histo1D(y_nbins, rap_beam - rap_max, rap_beam - rap_min);
+ for (size_t i = 0; i < 12; ++i) book(_tmphistos[i], "TMP/"+to_str(i), y_nbins, rap_min, rap_max);
+ for (size_t i = 12; i < 15; ++i) book(_tmphistos[i], "TMP/"+to_str(i), refData(dsShift+5, 1, 1));
+ for (size_t i = 15; i < 18; ++i) book(_tmphistos[i], "TMP/"+to_str(i), y_nbins, rap_beam - rap_max, rap_beam - rap_min);
+
+ int dsId = dsShift + 1;
+ for (size_t j = 0; j < 3; ++j) {
+ book(s1, dsId, 1, j+1);
+ book(s2, dsId+1, 1, j+1);
+ }
+ dsId += 2;
+ for (size_t j = 3; j < 6; ++j) {
+ book(s3, dsId, 1, 1);
+ dsId += 1;
+ book(s4, dsId, 1, 1);
+ dsId += 1;
+ }
declare(UnstableFinalState(), "UFS");
}
/// Perform the per-event analysis
void analyze(const Event& event) {
const double weight = 1.0;
const UnstableFinalState& ufs = apply<UnstableFinalState>(event, "UFS");
double ancestor_lftsum = 0.0;
double y, pT;
int id;
int partIdx = -1;
foreach (const Particle& p, ufs.particles()) {
id = p.pid();
// continue if particle not a K0s nor (anti-)Lambda
if ( (id == 310) || (id == -310) ) {
partIdx = 2;
} else if ( id == 3122 ) {
partIdx = 1;
} else if ( id == -3122 ) {
partIdx = 0;
} else {
continue;
}
ancestor_lftsum = getMotherLifeTimeSum(p);
// Lifetime cut: ctau sum of all particle ancestors < 10^-9 m according to the paper (see eq. 5)
const double MAX_CTAU = 1.0E-9; // [m]
if ( (ancestor_lftsum < 0.0) || (ancestor_lftsum > MAX_CTAU) ) continue;
const FourMomentum& qmom = p.momentum();
y = log((qmom.E() + qmom.pz())/(qmom.E() - qmom.pz()))/2.;
// skip this particle if it has too high or too low rapidity (extremely rare cases when E = +- pz)
if ( std::isnan(y) || std::isinf(y) ) continue;
y = fabs(y);
if (!inRange(y, rap_min, rap_max)) continue;
pT = sqrt((qmom.px() * qmom.px()) + (qmom.py() * qmom.py()));
if (!inRange(pT, pt_min, pt3_edge)) continue;
// Filling corresponding temporary histograms for pT intervals
- if (inRange(pT, pt_min, pt1_edge)) _tmphistos[partIdx*3].fill(y, weight);
- if (inRange(pT, pt1_edge, pt2_edge)) _tmphistos[partIdx*3+1].fill(y, weight);
- if (inRange(pT, pt2_edge, pt3_edge)) _tmphistos[partIdx*3+2].fill(y, weight);
+ if (inRange(pT, pt_min, pt1_edge)) _tmphistos[partIdx*3]->fill(y, weight);
+ if (inRange(pT, pt1_edge, pt2_edge)) _tmphistos[partIdx*3+1]->fill(y, weight);
+ if (inRange(pT, pt2_edge, pt3_edge)) _tmphistos[partIdx*3+2]->fill(y, weight);
// Fill histo in rapidity for whole pT interval
- _tmphistos[partIdx+9].fill(y, weight);
+ _tmphistos[partIdx+9]->fill(y, weight);
// Fill histo in pT for whole rapidity interval
- _tmphistos[partIdx+12].fill(pT, weight);
+ _tmphistos[partIdx+12]->fill(pT, weight);
// Fill histo in rapidity loss for whole pT interval
- _tmphistos[partIdx+15].fill(rap_beam - y, weight);
+ _tmphistos[partIdx+15]->fill(rap_beam - y, weight);
}
}
// Generate the ratio histograms
void finalize() {
int dsId = dsShift + 1;
for (size_t j = 0; j < 3; ++j) {
/// @todo Compactify to two one-liners
- Scatter2DPtr s1; book(s1, dsId, 1, j+1);
divide(_tmphistos[j], _tmphistos[3+j], s1);
- Scatter2DPtr s2; book(s2, dsId+1, 1, j+1);
divide(_tmphistos[j], _tmphistos[6+j], s2);
}
dsId += 2;
for (size_t j = 3; j < 6; ++j) {
/// @todo Compactify to two one-liners
- Scatter2DPtr s1; book(s1, dsId, 1, 1);
- divide(_tmphistos[3*j], _tmphistos[3*j+1], s1);
+ divide(_tmphistos[3*j], _tmphistos[3*j+1], s3);
dsId += 1;
- Scatter2DPtr s2; book(s2, dsId, 1, 1);
- divide(_tmphistos[3*j], _tmphistos[3*j+2], s2);
+ divide(_tmphistos[3*j], _tmphistos[3*j+2], s4);
dsId += 1;
}
}
//@}
private:
// Get particle lifetime from hardcoded data
double getLifeTime(int pid) {
double lft = -1.0;
if (pid < 0) pid = - pid;
// Correct Pythia6 PIDs for f0(980), f0(1370) mesons
if (pid == 10331) pid = 30221;
if (pid == 10221) pid = 9010221;
map<int, double>::iterator pPartLft = partLftMap.find(pid);
// search stable particle list
if (pPartLft == partLftMap.end()) {
if (pid <= 100) return 0.0;
for (size_t i=0; i < sizeof(stablePDGIds)/sizeof(unsigned int); i++) {
if (pid == stablePDGIds[i]) { lft = 0.0; break; }
}
} else {
lft = (*pPartLft).second;
}
if (lft < 0.0 && PID::isHadron(pid)) {
MSG_ERROR("Could not determine lifetime for particle with PID " << pid
<< "... This V^0 will be considered unprompt!");
}
return lft;
}
// Data members like post-cuts event weight counters go here
const double getMotherLifeTimeSum(const Particle& p) {
if (p.genParticle() == NULL) return -1.;
double lftSum = 0.;
double plft = 0.;
const GenParticle* part = p.genParticle();
const GenVertex* ivtx = part->production_vertex();
while (ivtx) {
if (ivtx->particles_in_size() < 1) { lftSum = -1.; break; };
const GenVertex::particles_in_const_iterator iPart_invtx = ivtx->particles_in_const_begin();
part = (*iPart_invtx);
if ( !(part) ) { lftSum = -1.; break; };
ivtx = part->production_vertex();
if ( (part->pdg_id() == 2212) || !(ivtx) ) break; //reached beam
plft = getLifeTime(part->pdg_id());
if (plft < 0.) { lftSum = -1.; break; };
lftSum += plft;
};
return (lftSum * c_light);
}
/// @name Private variables
//@{
// The rapidity of the beam according to the selected beam energy
double rap_beam;
// The edges of the intervals of transverse momentum
double pt_min, pt1_edge, pt2_edge, pt3_edge;
// The limits of the rapidity window
double rap_min;
double rap_max;
// Indicates which set of histograms will be output to yoda file (according to beam energy)
int dsShift;
// Map between PDG id and particle lifetimes in seconds
std::map<int, double> partLftMap;
// Set of PDG Ids for stable particles (PDG Id <= 100 are considered stable)
static const int stablePDGIds[205];
//@}
/// @name Helper histograms
//@{
/// Histograms are defined in the following order: anti-Lambda, Lambda and K0s.
/// First 3 suites of 3 histograms correspond to each particle in bins of y for the 3 pT intervals. (9 histos)
/// Next 3 histograms contain the particles in y bins for the whole pT interval (3 histos)
/// Next 3 histograms contain the particles in y_loss bins for the whole pT interval (3 histos)
/// Last 3 histograms contain the particles in pT bins for the whole rapidity (y) interval (3 histos)
- YODA::Histo1D _tmphistos[18];
+ Histo1DPtr _tmphistos[18];
+ Scatter2DPtr s1,s2,s3,s4;
//@}
// Fill the PDG Id to Lifetime[seconds] map
// Data was extracted from LHCb Particle Table through LHCb::ParticlePropertySvc
bool fillMap(map<int, double>& m) {
m[6] = 4.707703E-25; m[11] = 1.E+16; m[12] = 1.E+16;
m[13] = 2.197019E-06; m[14] = 1.E+16; m[15] = 2.906E-13; m[16] = 1.E+16; m[22] = 1.E+16;
m[23] = 2.637914E-25; m[24] = 3.075758E-25; m[25] = 9.4E-26; m[35] = 9.4E-26;
m[36] = 9.4E-26; m[37] = 9.4E-26; m[84] = 3.335641E-13; m[85] = 1.290893E-12;
m[111] = 8.4E-17; m[113] = 4.405704E-24; m[115] = 6.151516E-24; m[117] = 4.088275E-24;
m[119] = 2.102914E-24; m[130] = 5.116E-08; m[150] = 1.525E-12; m[211] = 2.6033E-08;
m[213] = 4.405704E-24; m[215] = 6.151516E-24; m[217] = 4.088275E-24; m[219] = 2.102914E-24;
m[221] = 5.063171E-19; m[223] = 7.752794E-23; m[225] = 3.555982E-24; m[227] = 3.91793E-24;
m[229] = 2.777267E-24; m[310] = 8.953E-11; m[313] = 1.308573E-23; m[315] = 6.038644E-24;
m[317] = 4.139699E-24; m[319] = 3.324304E-24; m[321] = 1.238E-08; m[323] = 1.295693E-23;
m[325] = 6.682357E-24; m[327] = 4.139699E-24; m[329] = 3.324304E-24; m[331] = 3.210791E-21;
m[333] = 1.545099E-22; m[335] = 9.016605E-24; m[337] = 7.565657E-24; m[350] = 1.407125E-12;
m[411] = 1.04E-12; m[413] = 6.856377E-21; m[415] = 1.778952E-23; m[421] = 4.101E-13;
m[423] = 1.000003E-19; m[425] = 1.530726E-23; m[431] = 5.E-13; m[433] = 1.000003E-19;
m[435] = 3.291061E-23; m[441] = 2.465214E-23; m[443] = 7.062363E-21; m[445] = 3.242425E-22;
m[510] = 1.525E-12; m[511] = 1.525E-12; m[513] = 1.000019E-19; m[515] = 1.31E-23;
m[521] = 1.638E-12; m[523] = 1.000019E-19; m[525] = 1.31E-23; m[530] = 1.536875E-12;
m[531] = 1.472E-12; m[533] = 1.E-19; m[535] = 1.31E-23; m[541] = 4.5E-13;
m[553] = 1.218911E-20; m[1112] = 4.539394E-24; m[1114] = 5.578069E-24; m[1116] = 1.994582E-24;
m[1118] = 2.269697E-24; m[1212] = 4.539394E-24; m[1214] = 5.723584E-24; m[1216] = 1.994582E-24;
m[1218] = 1.316424E-24; m[2112] = 8.857E+02; m[2114] = 5.578069E-24; m[2116] = 4.388081E-24;
m[2118] = 2.269697E-24; m[2122] = 4.539394E-24; m[2124] = 5.723584E-24; m[2126] = 1.994582E-24;
m[2128] = 1.316424E-24; m[2212] = 1.E+16; m[2214] = 5.578069E-24; m[2216] = 4.388081E-24;
m[2218] = 2.269697E-24; m[2222] = 4.539394E-24; m[2224] = 5.578069E-24; m[2226] = 1.994582E-24;
m[2228] = 2.269697E-24; m[3112] = 1.479E-10; m[3114] = 1.670589E-23; m[3116] = 5.485102E-24;
m[3118] = 3.656734E-24; m[3122] = 2.631E-10; m[3124] = 4.219309E-23; m[3126] = 8.227653E-24;
m[3128] = 3.291061E-24; m[3212] = 7.4E-20; m[3214] = 1.828367E-23; m[3216] = 5.485102E-24;
m[3218] = 3.656734E-24; m[3222] = 8.018E-11; m[3224] = 1.838582E-23; m[3226] = 5.485102E-24;
m[3228] = 3.656734E-24; m[3312] = 1.639E-10; m[3314] = 6.648608E-23; m[3322] = 2.9E-10;
m[3324] = 7.233101E-23; m[3334] = 8.21E-11; m[4112] = 2.991874E-22; m[4114] = 4.088274E-23;
m[4122] = 2.E-13; m[4132] = 1.12E-13; m[4212] = 3.999999E-22; m[4214] = 3.291061E-22;
m[4222] = 2.951624E-22; m[4224] = 4.417531E-23; m[4232] = 4.42E-13; m[4332] = 6.9E-14;
m[4412] = 3.335641E-13; m[4422] = 3.335641E-13; m[4432] = 3.335641E-13; m[5112] = 1.E-19;
m[5122] = 1.38E-12; m[5132] = 1.42E-12; m[5142] = 1.290893E-12; m[5212] = 1.E-19;
m[5222] = 1.E-19; m[5232] = 1.42E-12; m[5242] = 1.290893E-12; m[5312] = 1.E-19;
m[5322] = 1.E-19; m[5332] = 1.55E-12; m[5342] = 1.290893E-12; m[5442] = 1.290893E-12;
m[5512] = 1.290893E-12; m[5522] = 1.290893E-12; m[5532] = 1.290893E-12; m[5542] = 1.290893E-12;
m[10111] = 2.48382E-24; m[10113] = 4.635297E-24; m[10115] = 2.54136E-24; m[10211] = 2.48382E-24;
m[10213] = 4.635297E-24; m[10215] = 2.54136E-24; m[10223] = 1.828367E-24; m[10225] = 3.636531E-24;
m[10311] = 2.437823E-24; m[10313] = 7.313469E-24; m[10315] = 3.538775E-24;
m[10321] = 2.437823E-24; m[10323] = 7.313469E-24; m[10325] = 3.538775E-24;
m[10331] = 4.804469E-24; m[10411] = 4.38E-24; m[10413] = 3.29E-23; m[10421] = 4.38E-24;
m[10423] = 3.22653E-23; m[10431] = 6.5821E-22; m[10433] = 6.5821E-22; m[10441] = 6.453061E-23;
m[10511] = 4.39E-24; m[10513] = 1.65E-23; m[10521] = 4.39E-24; m[10523] = 1.65E-23;
m[10531] = 4.39E-24; m[10533] = 1.65E-23; m[11114] = 2.194041E-24; m[11116] = 1.828367E-24;
m[11212] = 1.880606E-24; m[11216] = 1.828367E-24; m[12112] = 2.194041E-24;
m[12114] = 2.194041E-24; m[12116] = 5.063171E-24; m[12126] = 1.828367E-24;
m[12212] = 2.194041E-24; m[12214] = 2.194041E-24; m[12216] = 5.063171E-24;
m[12224] = 2.194041E-24; m[12226] = 1.828367E-24; m[13112] = 6.582122E-24; m[13114] = 1.09702E-23;
m[13116] = 5.485102E-24; m[13122] = 1.316424E-23; m[13124] = 1.09702E-23; m[13126] = 6.928549E-24;
m[13212] = 6.582122E-24; m[13214] = 1.09702E-23; m[13216] = 5.485102E-24; m[13222] = 6.582122E-24;
m[13224] = 1.09702E-23; m[13226] = 5.485102E-24; m[13314] = 2.742551E-23;
m[13324] = 2.742551E-23; m[14122] = 1.828367E-22; m[20022] = 1.E+16; m[20113] = 1.567172E-24;
m[20213] = 1.567172E-24; m[20223] = 2.708692E-23; m[20313] = 3.782829E-24;
m[20315] = 2.384827E-24; m[20323] = 3.782829E-24; m[20325] = 2.384827E-24;
m[20333] = 1.198929E-23; m[20413] = 2.63E-24; m[20423] = 2.63E-24; m[20433] = 6.5821E-22;
m[20443] = 7.395643E-22; m[20513] = 2.63E-24; m[20523] = 2.63E-24; m[20533] = 2.63E-24;
m[21112] = 2.632849E-24; m[21114] = 3.291061E-24; m[21212] = 2.632849E-24;
m[21214] = 6.582122E-24; m[22112] = 4.388081E-24; m[22114] = 3.291061E-24;
m[22122] = 2.632849E-24; m[22124] = 6.582122E-24; m[22212] = 4.388081E-24;
m[22214] = 3.291061E-24; m[22222] = 2.632849E-24; m[22224] = 3.291061E-24;
m[23112] = 7.313469E-24; m[23114] = 2.991874E-24; m[23122] = 4.388081E-24;
m[23124] = 6.582122E-24; m[23126] = 3.291061E-24; m[23212] = 7.313469E-24;
m[23214] = 2.991874E-24; m[23222] = 7.313469E-24; m[23224] = 2.991874E-24;
m[30113] = 2.632849E-24; m[30213] = 2.632849E-24; m[30221] = 1.880606E-24;
m[30223] = 2.089563E-24; m[30313] = 2.056913E-24; m[30323] = 2.056913E-24;
m[30443] = 2.419898E-23; m[31114] = 1.880606E-24; m[31214] = 3.291061E-24;
m[32112] = 3.989164E-24; m[32114] = 1.880606E-24; m[32124] = 3.291061E-24;
m[32212] = 3.989164E-24; m[32214] = 1.880606E-24; m[32224] = 1.880606E-24;
m[33122] = 1.880606E-23; m[42112] = 6.582122E-24; m[42212] = 6.582122E-24;
m[43122] = 2.194041E-24; m[53122] = 4.388081E-24; m[100111] = 1.645531E-24;
m[100113] = 1.64553E-24; m[100211] = 1.645531E-24; m[100213] = 1.64553E-24;
m[100221] = 1.196749E-23; m[100223] = 3.061452E-24; m[100313] = 2.837122E-24;
m[100323] = 2.837122E-24; m[100331] = 4.459432E-25; m[100333] = 4.388081E-24;
m[100441] = 4.701516E-23; m[100443] = 2.076379E-21; m[100553] = 2.056913E-20;
m[200553] = 3.242425E-20; m[300553] = 3.210791E-23; m[9000111] = 8.776163E-24;
m[9000211] = 8.776163E-24; m[9000443] = 8.227652E-24; m[9000553] = 5.983747E-24;
m[9010111] = 3.164482E-24; m[9010211] = 3.164482E-24; m[9010221] = 9.403031E-24;
m[9010443] = 8.438618E-24; m[9010553] = 8.3318E-24; m[9020443] = 1.061633E-23;
m[9030221] = 6.038644E-24; m[9042413] = 2.07634E-21; m[9050225] = 1.394517E-24;
m[9060225] = 3.291061E-24; m[9080225] = 4.388081E-24; m[9090225] = 2.056913E-24;
m[9910445] = 2.07634E-21; m[9920443] = 2.07634E-21;
return true;
}
};
const int LHCB_2011_I917009::stablePDGIds[205] = {
311, 543, 545, 551, 555, 557, 1103, 2101, 2103, 2203, 3101, 3103, 3201, 3203, 3303,
4101, 4103, 4124, 4201, 4203, 4301, 4303, 4312, 4314, 4322, 4324, 4334, 4403, 4414,
4424, 4434, 4444, 5101, 5103, 5114, 5201, 5203, 5214, 5224, 5301, 5303, 5314, 5324,
5334, 5401, 5403, 5412, 5414, 5422, 5424, 5432, 5434, 5444, 5503, 5514, 5524, 5534,
5544, 5554, 10022, 10333, 10335, 10443, 10541, 10543, 10551, 10553, 10555, 11112,
12118, 12122, 12218, 12222, 13316, 13326, 20543, 20553, 20555, 23314, 23324, 30343,
30353, 30363, 30553, 33314, 33324, 41214, 42124, 52114, 52214, 100311, 100315, 100321,
100325, 100411, 100413, 100421, 100423, 100551, 100555, 100557, 110551, 110553, 110555,
120553, 120555, 130553, 200551, 200555, 210551, 210553, 220553, 1000001, 1000002,
1000003, 1000004, 1000005, 1000006, 1000011, 1000012, 1000013, 1000014, 1000015,
1000016, 1000021, 1000022, 1000023, 1000024, 1000025, 1000035, 1000037, 1000039,
2000001, 2000002, 2000003, 2000004, 2000005, 2000006, 2000011, 2000012, 2000013,
2000014, 2000015, 2000016, 3000111, 3000113, 3000211, 3000213, 3000221, 3000223,
3000331, 3100021, 3100111, 3100113, 3200111, 3200113, 3300113, 3400113, 4000001,
4000002, 4000011, 4000012, 5000039, 9000221, 9900012, 9900014, 9900016, 9900023,
9900024, 9900041, 9900042 };
// Hook for the plugin system
DECLARE_RIVET_PLUGIN(LHCB_2011_I917009);
}
+
diff --git a/src/Core/Analysis.cc b/src/Core/Analysis.cc
--- a/src/Core/Analysis.cc
+++ b/src/Core/Analysis.cc
@@ -1,862 +1,862 @@
// -*- C++ -*-
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Analysis.hh"
#include "Rivet/AnalysisHandler.hh"
#include "Rivet/AnalysisInfo.hh"
#include "Rivet/Tools/BeamConstraint.hh"
-#include "DummyConfig.hh"
-#ifdef HAVE_EXECINFO_H
-#include <execinfo.h>
-#endif
+// #include "DummyConfig.hh"
+// #ifdef HAVE_EXECINFO_H
+// #include <execinfo.h>
+// #endif
namespace Rivet {
Analysis::Analysis(const string& name)
: _crossSection(-1.0),
_gotCrossSection(false),
_analysishandler(NULL)
{
ProjectionApplier::_allowProjReg = false;
_defaultname = name;
unique_ptr<AnalysisInfo> ai = AnalysisInfo::make(name);
assert(ai);
_info = move(ai);
assert(_info);
}
double Analysis::sqrtS() const {
return handler().sqrtS();
}
const ParticlePair& Analysis::beams() const {
return handler().beams();
}
const PdgIdPair Analysis::beamIds() const {
return handler().beamIds();
}
const string Analysis::histoDir() const {
/// @todo Cache in a member variable
string _histoDir;
if (_histoDir.empty()) {
_histoDir = "/" + name();
if (handler().runName().length() > 0) {
_histoDir = "/" + handler().runName() + _histoDir;
}
replace_all(_histoDir, "//", "/"); //< iterates until none
}
return _histoDir;
}
const string Analysis::histoPath(const string& hname) const {
const string path = histoDir() + "/" + hname;
return path;
}
const string Analysis::histoPath(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
return histoDir() + "/" + makeAxisCode(datasetId, xAxisId, yAxisId);
}
const string Analysis::makeAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
stringstream axisCode;
axisCode << "d";
if (datasetId < 10) axisCode << 0;
axisCode << datasetId;
axisCode << "-x";
if (xAxisId < 10) axisCode << 0;
axisCode << xAxisId;
axisCode << "-y";
if (yAxisId < 10) axisCode << 0;
axisCode << yAxisId;
return axisCode.str();
}
Log& Analysis::getLog() const {
string logname = "Rivet.Analysis." + name();
return Log::getLog(logname);
}
size_t Analysis::numEvents() const {
return handler().numEvents();
}
double Analysis::sumOfWeights() const {
return handler().sumOfWeights();
}
///////////////////////////////////////////
bool Analysis::isCompatible(const ParticlePair& beams) const {
return isCompatible(beams.first.pid(), beams.second.pid(),
beams.first.energy(), beams.second.energy());
}
bool Analysis::isCompatible(PdgId beam1, PdgId beam2, double e1, double e2) const {
PdgIdPair beams(beam1, beam2);
pair<double,double> energies(e1, e2);
return isCompatible(beams, energies);
}
bool Analysis::isCompatible(const PdgIdPair& beams, const pair<double,double>& energies) const {
// First check the beam IDs
bool beamIdsOk = false;
foreach (const PdgIdPair& bp, requiredBeams()) {
if (compatible(beams, bp)) {
beamIdsOk = true;
break;
}
}
if (!beamIdsOk) return false;
// Next check that the energies are compatible (within 1% or 1 GeV, whichever is larger, for a bit of UI forgiveness)
/// @todo Use some sort of standard ordering to improve comparisons, esp. when the two beams are different particles
bool beamEnergiesOk = requiredEnergies().size() > 0 ? false : true;
typedef pair<double,double> DoublePair;
foreach (const DoublePair& ep, requiredEnergies()) {
if ((fuzzyEquals(ep.first, energies.first, 0.01) && fuzzyEquals(ep.second, energies.second, 0.01)) ||
(fuzzyEquals(ep.first, energies.second, 0.01) && fuzzyEquals(ep.second, energies.first, 0.01)) ||
(abs(ep.first - energies.first) < 1*GeV && abs(ep.second - energies.second) < 1*GeV) ||
(abs(ep.first - energies.second) < 1*GeV && abs(ep.second - energies.first) < 1*GeV)) {
beamEnergiesOk = true;
break;
}
}
return beamEnergiesOk;
/// @todo Need to also check internal consistency of the analysis'
/// beam requirements with those of the projections it uses.
}
///////////////////////////////////////////
Analysis& Analysis::setCrossSection(double xs) {
_crossSection = xs;
_gotCrossSection = true;
return *this;
}
double Analysis::crossSection() const {
if (!_gotCrossSection || std::isnan(_crossSection)) {
string errMsg = "You did not set the cross section for the analysis " + name();
throw Error(errMsg);
}
return _crossSection;
}
double Analysis::crossSectionPerEvent() const {
return _crossSection/sumOfWeights();
}
////////////////////////////////////////////////////////////
// Histogramming
void Analysis::_cacheRefData() const {
if (_refdata.empty()) {
MSG_TRACE("Getting refdata cache for paper " << name());
_refdata = getRefData(name());
}
}
CounterPtr & Analysis::book(CounterPtr & ctr,
const string& cname,
const string& title) {
const string path = histoPath(cname);
ctr = CounterPtr(handler().numWeights(), Counter(path, title));
addAnalysisObject(ctr);
MSG_TRACE("Made counter " << cname << " for " << name());
return ctr;
}
CounterPtr & Analysis::book(CounterPtr & ctr, unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
const string& title) {
// const string& xtitle,
// const string& ytitle) {
const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
return book(ctr, axisCode, title);
}
Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname,
size_t nbins, double lower, double upper,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
Histo1D hist = Histo1D(nbins, lower, upper, path, title);
hist.setAnnotation("XLabel", xtitle);
hist.setAnnotation("YLabel", ytitle);
histo = Histo1DPtr(handler().numWeights(), hist);
addAnalysisObject(histo);
MSG_TRACE("Made histogram " << hname << " for " << name());
return histo;
// Histo1DPtr hist;
// try { // try to bind to pre-existing
// // AnalysisObjectPtr ao = getAnalysisObject(path);
// // hist = dynamic_pointer_cast<Histo1D>(ao);
// hist = getHisto1D(hname);
// /// @todo Test that cast worked
// /// @todo Also test that binning is as expected?
// MSG_TRACE("Bound pre-existing histogram " << hname << " for " << name());
// } catch (...) { // binding failed; make it from scratch
// hist = make_shared<Histo1D>(nbins, lower, upper, histoPath(hname), title);
// addAnalysisObject(hist);
// MSG_TRACE("Made histogram " << hname << " for " << name());
// }
}
Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname,
const initializer_list<double>& binedges,
const string& title,
const string& xtitle,
const string& ytitle) {
return book(histo, hname, vector<double>{binedges}, title, xtitle, ytitle);
}
Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname,
const vector<double>& binedges,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
// Histo1DPtr hist;
// try { // try to bind to pre-existing
// // AnalysisObjectPtr ao = getAnalysisObject(path);
// // hist = dynamic_pointer_cast<Histo1D>(ao);
// hist = getHisto1D(hname);
// /// @todo Test that cast worked
// /// @todo Also test that binning is as expected?
// MSG_TRACE("Bound pre-existing histogram " << hname << " for " << name());
// } catch (...) { // binding failed; make it from scratch
// hist = make_shared<Histo1D>(binedges, histoPath(hname), title);
// addAnalysisObject(hist);
// MSG_TRACE("Made histogram " << hname << " for " << name());
// }
Histo1D hist = Histo1D(binedges, path, title);
hist.setAnnotation("XLabel", xtitle);
hist.setAnnotation("YLabel", ytitle);
histo = Histo1DPtr(handler().numWeights(), hist);
addAnalysisObject(histo);
MSG_TRACE("Made histogram " << hname << " for " << name());
return histo;
}
Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname,
const Scatter2D& refscatter,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
// Histo1DPtr hist;
// try { // try to bind to pre-existing
// // AnalysisObjectPtr ao = getAnalysisObject(path);
// // hist = dynamic_pointer_cast<Histo1D>(ao);
// hist = getHisto1D(hname);
// /// @todo Test that cast worked
// /// @todo Also test that binning is as expected?
// MSG_TRACE("Bound pre-existing histogram " << hname << " for " << name());
// } catch (...) { // binding failed; make it from scratch
// hist = make_shared<Histo1D>(refscatter, histoPath(hname));
// if (hist->hasAnnotation("IsRef")) hist->rmAnnotation("IsRef");
// addAnalysisObject(hist);
// MSG_TRACE("Made histogram " << hname << " for " << name());
// }
Histo1D hist = Histo1D(refscatter, path);
hist.setTitle(title);
hist.setAnnotation("XLabel", xtitle);
hist.setAnnotation("YLabel", ytitle);
histo = Histo1DPtr(handler().numWeights(), hist);
addAnalysisObject(histo);
MSG_TRACE("Made histogram " << hname << " for " << name());
return histo;
}
Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname,
const string& title,
const string& xtitle,
const string& ytitle) {
const Scatter2D& refdata = refData(hname);
return book(histo, hname, refdata, title, xtitle, ytitle);
}
Histo1DPtr & Analysis::book(Histo1DPtr & histo, unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
const string& title,
const string& xtitle,
const string& ytitle) {
const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
return book(histo, axisCode, title, xtitle, ytitle);
}
/// @todo Add booking methods which take a path, titles and *a reference Scatter from which to book*
/////////////////
Histo2DPtr & Analysis::book(Histo2DPtr & h2d,const string& hname,
size_t nxbins, double xlower, double xupper,
size_t nybins, double ylower, double yupper,
const string& title,
const string& xtitle,
const string& ytitle,
const string& ztitle)
{
const string path = histoPath(hname);
Histo2D hist(nxbins, xlower, xupper, nybins, ylower, yupper, path, title);
hist.setAnnotation("XLabel", xtitle);
hist.setAnnotation("YLabel", ytitle);
hist.setAnnotation("ZLabel", ztitle);
h2d = Histo2DPtr(handler().numWeights(), hist);
addAnalysisObject(h2d);
MSG_TRACE("Made 2D histogram " << hname << " for " << name());
return h2d;
}
Histo2DPtr & Analysis::book(Histo2DPtr & h2d,const string& hname,
const initializer_list<double>& xbinedges,
const initializer_list<double>& ybinedges,
const string& title,
const string& xtitle,
const string& ytitle,
const string& ztitle)
{
return book(h2d, hname, vector<double>{xbinedges}, vector<double>{ybinedges}, title, xtitle, ytitle, ztitle);
}
Histo2DPtr & Analysis::book(Histo2DPtr & h2d,const string& hname,
const vector<double>& xbinedges,
const vector<double>& ybinedges,
const string& title,
const string& xtitle,
const string& ytitle,
const string& ztitle)
{
const string path = histoPath(hname);
Histo2D hist(xbinedges, ybinedges, path, title);
hist.setAnnotation("XLabel", xtitle);
hist.setAnnotation("YLabel", ytitle);
hist.setAnnotation("ZLabel", ztitle);
h2d = Histo2DPtr(handler().numWeights(), hist);
addAnalysisObject(h2d);
MSG_TRACE("Made 2D histogram " << hname << " for " << name());
return h2d;
}
Profile1DPtr & Analysis::book(Profile1DPtr & p1d,const string& hname,
size_t nbins, double lower, double upper,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
Profile1D prof(nbins, lower, upper, path, title);
prof.setAnnotation("XLabel", xtitle);
prof.setAnnotation("YLabel", ytitle);
p1d = Profile1DPtr(handler().numWeights(), prof);
addAnalysisObject(p1d);
MSG_TRACE("Made profile histogram " << hname << " for " << name());
return p1d;
}
Profile1DPtr & Analysis::book(Profile1DPtr & p1d,const string& hname,
const initializer_list<double>& binedges,
const string& title,
const string& xtitle,
const string& ytitle) {
return book(p1d, hname, vector<double>{binedges}, title, xtitle, ytitle);
}
Profile1DPtr & Analysis::book(Profile1DPtr & p1d,const string& hname,
const vector<double>& binedges,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
Profile1D prof(binedges, path, title);
prof.setAnnotation("XLabel", xtitle);
prof.setAnnotation("YLabel", ytitle);
p1d = Profile1DPtr(handler().numWeights(), prof);
addAnalysisObject(p1d);
MSG_TRACE("Made profile histogram " << hname << " for " << name());
return p1d;
}
Profile1DPtr & Analysis::book(Profile1DPtr & p1d,const string& hname,
const Scatter2D& refscatter,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
Profile1D prof(refscatter, path);
prof.setTitle(title);
prof.setAnnotation("XLabel", xtitle);
prof.setAnnotation("YLabel", ytitle);
p1d = Profile1DPtr(handler().numWeights(), prof);
addAnalysisObject(p1d);
MSG_TRACE("Made profile histogram " << hname << " for " << name());
return p1d;
// if (prof.hasAnnotation("IsRef")) prof.rmAnnotation("IsRef");
}
Profile1DPtr & Analysis::book(Profile1DPtr & p1d,const string& hname,
const string& title,
const string& xtitle,
const string& ytitle) {
const Scatter2D& refdata = refData(hname);
book(p1d, hname, refdata, title, xtitle, ytitle);
return p1d;
}
Profile1DPtr & Analysis::book(Profile1DPtr & p1d,unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
const string& title,
const string& xtitle,
const string& ytitle) {
const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
return book(p1d, axisCode, title, xtitle, ytitle);
}
Profile2DPtr & Analysis::book(Profile2DPtr & p2d, const string& hname,
size_t nxbins, double xlower, double xupper,
size_t nybins, double ylower, double yupper,
const string& title,
const string& xtitle,
const string& ytitle,
const string& ztitle)
{
const string path = histoPath(hname);
Profile2D prof(nxbins, xlower, xupper, nybins, ylower, yupper, path, title);
prof.setAnnotation("XLabel", xtitle);
prof.setAnnotation("YLabel", ytitle);
prof.setAnnotation("ZLabel", ztitle);
p2d = Profile2DPtr(handler().numWeights(), prof);
addAnalysisObject(p2d);
MSG_TRACE("Made 2D profile histogram " << hname << " for " << name());
return p2d;
}
Profile2DPtr & Analysis::book(Profile2DPtr & p2d, const string& hname,
const initializer_list<double>& xbinedges,
const initializer_list<double>& ybinedges,
const string& title,
const string& xtitle,
const string& ytitle,
const string& ztitle)
{
return book(p2d, hname, vector<double>{xbinedges}, vector<double>{ybinedges}, title, xtitle, ytitle, ztitle);
}
Profile2DPtr & Analysis::book(Profile2DPtr & p2d, const string& hname,
const vector<double>& xbinedges,
const vector<double>& ybinedges,
const string& title,
const string& xtitle,
const string& ytitle,
const string& ztitle)
{
const string path = histoPath(hname);
Profile2D prof(xbinedges, ybinedges, path, title);
prof.setAnnotation("XLabel", xtitle);
prof.setAnnotation("YLabel", ytitle);
prof.setAnnotation("ZLabel", ztitle);
p2d = Profile2DPtr(handler().numWeights(), prof);
addAnalysisObject(p2d);
MSG_TRACE("Made 2D profile histogram " << hname << " for " << name());
return p2d;
}
Scatter2DPtr & Analysis::book(Scatter2DPtr & s2d, unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
bool copy_pts,
const string& title,
const string& xtitle,
const string& ytitle) {
const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
return book(s2d, axisCode, copy_pts, title, xtitle, ytitle);
}
Scatter2DPtr & Analysis::book(Scatter2DPtr & s2d, const string& hname,
bool copy_pts,
const string& title,
const string& xtitle,
const string& ytitle) {
Scatter2D scat;
const string path = histoPath(hname);
if (copy_pts) {
const Scatter2D& refdata = refData(hname);
scat = Scatter2D(refdata, path);
for (Point2D& p : scat.points()) p.setY(0, 0);
} else {
scat = Scatter2D(path);
}
scat.setTitle(title);
scat.setAnnotation("XLabel", xtitle);
scat.setAnnotation("YLabel", ytitle);
s2d = Scatter2DPtr(handler().numWeights(), scat);
addAnalysisObject(s2d);
MSG_TRACE("Made scatter " << hname << " for " << name());
// if (scat.hasAnnotation("IsRef")) scat.rmAnnotation("IsRef");
return s2d;
}
Scatter2DPtr & Analysis::book(Scatter2DPtr & s2d, const string& hname,
size_t npts, double lower, double upper,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
Scatter2D scat;
const double binwidth = (upper-lower)/npts;
for (size_t pt = 0; pt < npts; ++pt) {
const double bincentre = lower + (pt + 0.5) * binwidth;
scat.addPoint(bincentre, 0, binwidth/2.0, 0);
}
scat.setTitle(title);
scat.setAnnotation("XLabel", xtitle);
scat.setAnnotation("YLabel", ytitle);
s2d = Scatter2DPtr(handler().numWeights(), scat);
addAnalysisObject(s2d);
MSG_TRACE("Made scatter " << hname << " for " << name());
return s2d;
}
Scatter2DPtr & Analysis::book(Scatter2DPtr & s2d, const string& hname,
const vector<double>& binedges,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
Scatter2D scat;
for (size_t pt = 0; pt < binedges.size()-1; ++pt) {
const double bincentre = (binedges[pt] + binedges[pt+1]) / 2.0;
const double binwidth = binedges[pt+1] - binedges[pt];
scat.addPoint(bincentre, 0, binwidth/2.0, 0);
}
scat.setTitle(title);
scat.setAnnotation("XLabel", xtitle);
scat.setAnnotation("YLabel", ytitle);
s2d = Scatter2DPtr(handler().numWeights(), scat);
addAnalysisObject(s2d);
MSG_TRACE("Made scatter " << hname << " for " << name());
return s2d;
}
void Analysis::divide(CounterPtr c1, CounterPtr c2, Scatter1DPtr s) const {
const string path = s->path();
*s = *c1 / *c2;
s->setPath(path);
}
void Analysis::divide(const Counter& c1, const Counter& c2, Scatter1DPtr s) const {
const string path = s->path();
*s = c1 / c2;
s->setPath(path);
}
void Analysis::divide(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const {
const string path = s->path();
*s = *h1 / *h2;
s->setPath(path);
}
void Analysis::divide(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const {
const string path = s->path();
*s = h1 / h2;
s->setPath(path);
}
void Analysis::divide(Profile1DPtr p1, Profile1DPtr p2, Scatter2DPtr s) const {
const string path = s->path();
*s = *p1 / *p2;
s->setPath(path);
}
void Analysis::divide(const Profile1D& p1, const Profile1D& p2, Scatter2DPtr s) const {
const string path = s->path();
*s = p1 / p2;
s->setPath(path);
}
void Analysis::divide(Histo2DPtr h1, Histo2DPtr h2, Scatter3DPtr s) const {
const string path = s->path();
*s = *h1 / *h2;
s->setPath(path);
}
void Analysis::divide(const Histo2D& h1, const Histo2D& h2, Scatter3DPtr s) const {
const string path = s->path();
*s = h1 / h2;
s->setPath(path);
}
void Analysis::divide(Profile2DPtr p1, Profile2DPtr p2, Scatter3DPtr s) const {
const string path = s->path();
*s = *p1 / *p2;
s->setPath(path);
}
void Analysis::divide(const Profile2D& p1, const Profile2D& p2, Scatter3DPtr s) const {
const string path = s->path();
*s = p1 / p2;
s->setPath(path);
}
/// @todo Counter and Histo2D efficiencies and asymms
void Analysis::efficiency(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const {
const string path = s->path();
*s = YODA::efficiency(*h1, *h2);
s->setPath(path);
}
void Analysis::efficiency(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const {
const string path = s->path();
*s = YODA::efficiency(h1, h2);
s->setPath(path);
}
void Analysis::asymm(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const {
const string path = s->path();
*s = YODA::asymm(*h1, *h2);
s->setPath(path);
}
void Analysis::asymm(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const {
const string path = s->path();
*s = YODA::asymm(h1, h2);
s->setPath(path);
}
void Analysis::scale(CounterPtr cnt, double factor) {
if (!cnt) {
MSG_WARNING("Failed to scale counter=NULL in analysis " << name() << " (scale=" << factor << ")");
return;
}
if (std::isnan(factor) || std::isinf(factor)) {
MSG_WARNING("Failed to scale counter=" << cnt->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")");
factor = 0;
}
MSG_TRACE("Scaling counter " << cnt->path() << " by factor " << factor);
try {
cnt->scaleW(factor);
} catch (YODA::Exception& we) {
MSG_WARNING("Could not scale counter " << cnt->path());
return;
}
}
void Analysis::normalize(Histo1DPtr histo, double norm, bool includeoverflows) {
if (!histo) {
MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << norm << ")");
return;
}
MSG_TRACE("Normalizing histo " << histo->path() << " to " << norm);
try {
histo->normalize(norm, includeoverflows);
} catch (YODA::Exception& we) {
MSG_WARNING("Could not normalize histo " << histo->path());
return;
}
}
void Analysis::scale(Histo1DPtr histo, double factor) {
if (!histo) {
MSG_WARNING("Failed to scale histo=NULL in analysis " << name() << " (scale=" << factor << ")");
return;
}
if (std::isnan(factor) || std::isinf(factor)) {
MSG_WARNING("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")");
factor = 0;
}
MSG_TRACE("Scaling histo " << histo->path() << " by factor " << factor);
try {
histo->scaleW(factor);
} catch (YODA::Exception& we) {
MSG_WARNING("Could not scale histo " << histo->path());
return;
}
}
void Analysis::normalize(Histo2DPtr histo, double norm, bool includeoverflows) {
if (!histo) {
MSG_ERROR("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << norm << ")");
return;
}
MSG_TRACE("Normalizing histo " << histo->path() << " to " << norm);
try {
histo->normalize(norm, includeoverflows);
} catch (YODA::Exception& we) {
MSG_WARNING("Could not normalize histo " << histo->path());
return;
}
}
void Analysis::scale(Histo2DPtr histo, double factor) {
if (!histo) {
MSG_ERROR("Failed to scale histo=NULL in analysis " << name() << " (scale=" << factor << ")");
return;
}
if (std::isnan(factor) || std::isinf(factor)) {
MSG_ERROR("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")");
factor = 0;
}
MSG_TRACE("Scaling histo " << histo->path() << " by factor " << factor);
try {
histo->scaleW(factor);
} catch (YODA::Exception& we) {
MSG_WARNING("Could not scale histo " << histo->path());
return;
}
}
void Analysis::integrate(Histo1DPtr h, Scatter2DPtr s) const {
// preserve the path info
const string path = s->path();
*s = toIntegralHisto(*h);
s->setPath(path);
}
void Analysis::integrate(const Histo1D& h, Scatter2DPtr s) const {
// preserve the path info
const string path = s->path();
*s = toIntegralHisto(h);
s->setPath(path);
}
}
/// @todo 2D versions of integrate... defined how, exactly?!?
//////////////////////////////////
namespace {
void errormsg(std::string name) {
-#ifdef HAVE_BACKTRACE
- void * buffer[4];
- backtrace(buffer, 4);
- backtrace_symbols_fd(buffer, 4 , 1);
-#endif
+// #ifdef HAVE_BACKTRACE
+// void * buffer[4];
+// backtrace(buffer, 4);
+// backtrace_symbols_fd(buffer, 4 , 1);
+// #endif
std::cerr << name << ": Can't book objects outside of init().\n";
assert(false);
}
}
namespace Rivet {
void Analysis::addAnalysisObject(const MultiweightAOPtr & ao) {
if (handler().stage() == AnalysisHandler::Stage::INIT) {
_analysisobjects.push_back(ao);
ao.get()->blockDestructor(true);
}
else {
errormsg(name());
}
}
void Analysis::removeAnalysisObject(const string& path) {
for (auto it = _analysisobjects.begin();
it != _analysisobjects.end(); ++it) {
if ((*it)->path() == path) {
_analysisobjects.erase(it);
break;
}
}
}
void Analysis::removeAnalysisObject(const MultiweightAOPtr & ao) {
for (auto it = _analysisobjects.begin(); it != _analysisobjects.end(); ++it) {
if ((*it) == ao) {
_analysisobjects.erase(it);
break;
}
}
}
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:43 PM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805897
Default Alt Text
(76 KB)

Event Timeline