|
|
[Rivet-svn] r2185 - in trunk: data/plotinfo src/Analyses src/Projections
blackhole at projects.hepforge.org
blackhole at projects.hepforge.org
Sun Dec 13 16:39:36 GMT 2009
Author: buckley
Date: Sun Dec 13 16:39:36 2009
New Revision: 2185
Log:
A few improvements to the W validation analysis
Modified:
trunk/data/plotinfo/MC_WANALYSIS.plot
trunk/src/Analyses/MC_WANALYSIS.cc
trunk/src/Projections/InvMassFinalState.cc
trunk/src/Projections/VetoedFinalState.cc
trunk/src/Projections/WFinder.cc
Modified: trunk/data/plotinfo/MC_WANALYSIS.plot
==============================================================================
--- trunk/data/plotinfo/MC_WANALYSIS.plot Sun Dec 13 16:29:51 2009 (r2184)
+++ trunk/data/plotinfo/MC_WANALYSIS.plot Sun Dec 13 16:39:36 2009 (r2185)
@@ -7,21 +7,21 @@
# BEGIN PLOT /MC_WANALYSIS/w-n
Title=W candidate multiplicity
-XLabel=$n_\text{W}^\pm$
+XLabel=$n(\text{W}^\pm)$
YLabel=$\mathrm{d}N_\text{evt}/\mathrm{d}n_\text{W}$
LogY=0
# END PLOT
# BEGIN PLOT /MC_WANALYSIS/wplus-n
Title=\PWplus candidate multiplicity
-XLabel=$n_\text{W}^+$
+XLabel=$n(\text{W}^+)$
YLabel=$\mathrm{d}N_\text{evt}/\mathrm{d}n_\text{W}$
LogY=0
# END PLOT
# BEGIN PLOT /MC_WANALYSIS/wminus-n
Title=\PWminus candidate multiplicity
-XLabel=$n_\text{W}^-$
+XLabel=$n(\text{W}^-)$
YLabel=$\mathrm{d}N_\text{evt}/\mathrm{d}n_\text{W}$
LogY=0
# END PLOT
Modified: trunk/src/Analyses/MC_WANALYSIS.cc
==============================================================================
--- trunk/src/Analyses/MC_WANALYSIS.cc Sun Dec 13 16:29:51 2009 (r2184)
+++ trunk/src/Analyses/MC_WANALYSIS.cc Sun Dec 13 16:39:36 2009 (r2185)
@@ -16,6 +16,9 @@
MC_WANALYSIS() : Analysis("MC_WANALYSIS")
{
setNeedsCrossSection(true);
+ _sumw = 0;
+ _sumw_wplus = 0;
+ _sumw_wminus = 0;
}
@@ -24,15 +27,17 @@
//@{
void init() {
+ // Projections
const ChargedFinalState cfs(-2, 2, 200*MeV);
addProjection(cfs, "CFS");
-
- /// @todo Handle muon-decay Ws as well
- const WFinder wf(-2, 2, 10.0*GeV, ELECTRON, 60.0*GeV, 100.0*GeV, 0.2);
- addProjection(wf, "WF");
- FastJets fastjets(wf.remainingFinalState(), FastJets::KT, 1.0);
+ const WFinder wfe(-2, 2, 10.0*GeV, ELECTRON, 60.0*GeV, 100.0*GeV, 0.2);
+ addProjection(wfe, "WFe");
+ const WFinder wfmu(-2, 2, 10.0*GeV, MUON, 60.0*GeV, 100.0*GeV, 0.2);
+ addProjection(wfmu, "WFmu");
+ FastJets fastjets(wfe.remainingFinalState(), FastJets::KT, 0.5);
addProjection(fastjets, "Jets");
+ // Histos
_hist_chargemulti = bookHistogram1D("track-n", 50, -0.5, 999.5);
_hist_chargept = bookHistogram1D("track-pt", 20, 0, 20);
/// @todo Use profile plots instead:
@@ -60,16 +65,18 @@
_hist_jetcount = bookHistogram1D("jet-n", 6, -0.5, 5.5);
_hist_jetpt = bookHistogram1D("jet-pt", 50, 20, 100);
}
+
void analyze(const Event& event) {
- const WFinder& wf = applyProjection<WFinder>(event, "WF");
+ const WFinder& wf = applyProjection<WFinder>(event, "WFe");
if (wf.size() == 0) {
getLog() << Log::DEBUG << "No W candidates found: vetoing" << endl;
vetoEvent;
}
const double weight = event.weight();
-
+ _sumw += weight;
+
// Charged particles part
const FinalState& cfs = applyProjection<FinalState>(event, "CFS");
_hist_chargemulti->fill(cfs.particles().size(), weight);
@@ -84,7 +91,7 @@
_hist_chargemeanpt->fill(meanpt/GeV, weight);
rmspt = sqrt(rmspt / cfs.particles().size());
_hist_chargermspt->fill(rmspt/GeV, weight);
-
+
// W part
unsigned int n_wplus(0), n_wminus(0);
foreach (const Particle& wp, wf.particles()) {
@@ -99,12 +106,14 @@
_hist_wmass->fill(m/GeV, weight);
if (wp.pdgId() == WPLUSBOSON) {
n_wplus += 1;
+ _sumw_wplus += weight;
_hist_wpluspt->fill(pT/GeV, weight);
_hist_wpluseta->fill(eta, weight);
_hist_wplusphi->fill(phi, weight);
_hist_wplusmass->fill(m/GeV, weight);
} else if (wp.pdgId() == WMINUSBOSON) {
n_wminus += 1;
+ _sumw_wminus += weight;
_hist_wminuspt->fill(pT/GeV, weight);
_hist_wminuseta->fill(eta, weight);
_hist_wminusphi->fill(phi, weight);
@@ -121,7 +130,6 @@
// Jet part
const FastJets& fastjets = applyProjection<FastJets>(event, "Jets");
const Jets jets = fastjets.jetsByPt(10*GeV);
- cout << jets.size() << endl;
_hist_jetcount->fill(jets.size(), weight);
foreach (const Jet& j, jets) {
const double pT = j.momentum().pT();
@@ -133,9 +141,8 @@
void finalize() {
const double xsec_sumw = crossSectionPerEvent()/picobarn;
- /// @todo Actually "measure" separate W+ and W- xsecs
- const double xsec_sumw_plus = xsec_sumw/2.0;
- const double xsec_sumw_minus = xsec_sumw/2.0;
+ const double xsec_sumw_plus = _sumw_wplus/_sumw * xsec_sumw;
+ const double xsec_sumw_minus = _sumw_wminus/_sumw * xsec_sumw;
scale(_hist_chargemulti, xsec_sumw);
scale(_hist_chargept, xsec_sumw);
@@ -165,6 +172,13 @@
private:
+ /// @name Weight counters
+ //@{
+ double _sumw;
+ double _sumw_wplus;
+ double _sumw_wminus;
+ //@}
+
/// @name Histograms
//@{
AIDA::IHistogram1D *_hist_chargemulti;
Modified: trunk/src/Projections/InvMassFinalState.cc
==============================================================================
--- trunk/src/Projections/InvMassFinalState.cc Sun Dec 13 16:29:51 2009 (r2184)
+++ trunk/src/Projections/InvMassFinalState.cc Sun Dec 13 16:39:36 2009 (r2185)
@@ -35,7 +35,7 @@
if (fscmp != EQUIVALENT) return fscmp;
// Then compare the two as final states
- const InvMassFinalState & other = dynamic_cast <const InvMassFinalState&>(p);
+ const InvMassFinalState& other = dynamic_cast<const InvMassFinalState&>(p);
fscmp = FinalState::compare(other);
if (fscmp != EQUIVALENT) return fscmp;
@@ -56,8 +56,8 @@
void InvMassFinalState::project(const Event& e) {
- const FinalState& fs = applyProjection<FinalState>(e, "FS");
_theParticles.clear();
+ const FinalState& fs = applyProjection<FinalState>(e, "FS");
// Containers for the particles of type specified in the pair
vector<const Particle*> type1;
Modified: trunk/src/Projections/VetoedFinalState.cc
==============================================================================
--- trunk/src/Projections/VetoedFinalState.cc Sun Dec 13 16:29:51 2009 (r2184)
+++ trunk/src/Projections/VetoedFinalState.cc Sun Dec 13 16:39:36 2009 (r2185)
@@ -42,14 +42,14 @@
codes.push_back(code->first);
}
const string codestr = "{ " + join(codes) + " }";
- getLog() << Log::DEBUG << p.pdgId() << " vs. veto codes = "
+ getLog() << Log::TRACE << p.pdgId() << " vs. veto codes = "
<< codestr << " (" << codes.size() << ")" << endl;
}
const long pdgid = p.pdgId();
const double pt = p.momentum().pT();
VetoDetails::iterator iter = _vetoCodes.find(pdgid);
if (iter == _vetoCodes.end()) {
- getLog() << Log::DEBUG << "Storing with PDG code = " << pdgid << ", pT = " << pt << endl;
+ getLog() << Log::TRACE << "Storing with PDG code = " << pdgid << ", pT = " << pt << endl;
_theParticles.push_back(p);
} else {
// This particle code is listed as a possible veto... check pT.
@@ -60,14 +60,14 @@
if (ptrange.first < numeric_limits<double>::max()) rangess << ptrange.second;
rangess << " - ";
if (ptrange.second < numeric_limits<double>::max()) rangess << ptrange.second;
- getLog() << Log::DEBUG << "ID = " << pdgid << ", pT range = " << rangess.str();
+ getLog() << Log::TRACE << "ID = " << pdgid << ", pT range = " << rangess.str();
stringstream debugline;
debugline << "with PDG code = " << pdgid << " pT = " << p.momentum().pT();
if (pt < ptrange.first || pt > ptrange.second) {
- getLog() << Log::DEBUG << "Storing " << debugline.str() << endl;
+ getLog() << Log::TRACE << "Storing " << debugline.str() << endl;
_theParticles.push_back(p);
} else {
- getLog() << Log::DEBUG << "Vetoing " << debugline.str() << endl;
+ getLog() << Log::TRACE << "Vetoing " << debugline.str() << endl;
}
}
}
@@ -118,16 +118,15 @@
for (set<ParticleVector::iterator>::reverse_iterator p = toErase.rbegin(); p != toErase.rend(); ++p) {
_theParticles.erase(*p);
}
-
+
+ /// @todo Improve!
for (ParentVetos::const_iterator vIt = _parentVetoes.begin(); vIt != _parentVetoes.end(); ++vIt) {
for (ParticleVector::iterator p = _theParticles.begin(); p != _theParticles.end(); ++p) {
- GenVertex *startVtx=((*p).genParticle()).production_vertex();
+ GenVertex* startVtx = ((*p).genParticle()).production_vertex();
bool veto = false;
- GenParticle HepMCP = (*p).genParticle();
if (startVtx!=0) {
for (GenVertex::particle_iterator pIt = startVtx->particles_begin(HepMC::ancestors);
- pIt != startVtx->particles_end(HepMC::ancestors) && !veto; ++pIt) {
-
+ pIt != startVtx->particles_end(HepMC::ancestors) && !veto; ++pIt) {
if (*vIt == (*pIt)->pdg_id()) {
veto = true;
p = _theParticles.erase(p);
@@ -137,17 +136,18 @@
}
}
}
-
+
// Now veto on the FS
foreach (const string& ifs, _vetofsnames) {
const FinalState& vfs = applyProjection<FinalState>(e, ifs);
const ParticleVector& vfsp = vfs.particles();
- for (ParticleVector::iterator icheck = _theParticles.begin(); icheck != _theParticles.end(); ++icheck){
+ for (ParticleVector::iterator icheck = _theParticles.begin(); icheck != _theParticles.end(); ++icheck) {
if (!icheck->hasGenParticle()) continue;
bool found = false;
for (ParticleVector::const_iterator ipart = vfsp.begin(); ipart != vfsp.end(); ++ipart){
if (!ipart->hasGenParticle()) continue;
- getLog() << Log::DEBUG << "Comparing barcode " << icheck->genParticle().barcode() << " with veto particle " << ipart->genParticle().barcode() << endl;
+ getLog() << Log::TRACE << "Comparing barcode " << icheck->genParticle().barcode()
+ << " with veto particle " << ipart->genParticle().barcode() << endl;
if (ipart->genParticle().barcode() == icheck->genParticle().barcode()){
found = true;
break;
Modified: trunk/src/Projections/WFinder.cc
==============================================================================
--- trunk/src/Projections/WFinder.cc Sun Dec 13 16:29:51 2009 (r2184)
+++ trunk/src/Projections/WFinder.cc Sun Dec 13 16:39:36 2009 (r2185)
@@ -32,14 +32,15 @@
WFinder::WFinder(const std::vector<std::pair<double, double> >& etaRanges,
double pTmin,
PdgId pid,
- double m2_min, const double m2_max,
+ double m2_min, double m2_max,
double dRmax) {
_init(etaRanges, pTmin, pid, m2_min, m2_max, dRmax);
}
void WFinder::_init(const std::vector<std::pair<double, double> >& etaRanges,
- double pTmin, PdgId pid,
+ double pTmin,
+ PdgId pid,
double m2_min, double m2_max,
double dRmax) {
FinalState fs(etaRanges, pTmin);
@@ -54,14 +55,14 @@
{
setName("WFinder");
- addProjection(fs, "FS");
+ //addProjection(fs, "FS");
assert(abs(pid) == ELECTRON || abs(pid) == MUON || abs(pid) == TAU);
PdgId nu_pid = abs(pid) + 1;
assert(abs(nu_pid) == NU_E || abs(nu_pid) == NU_MU || abs(nu_pid) == NU_TAU);
- std::vector<std::pair<long, long> > l_nu_ids;
- l_nu_ids += std::make_pair(abs(pid), -abs(nu_pid));
- l_nu_ids += std::make_pair(-abs(pid), abs(nu_pid));
+ vector<pair<long, long> > l_nu_ids;
+ l_nu_ids += make_pair(abs(pid), -abs(nu_pid));
+ l_nu_ids += make_pair(-abs(pid), abs(nu_pid));
InvMassFinalState imfs(fs, l_nu_ids, m2_min, m2_max);
addProjection(imfs, "IMFS");
More information about the Rivet-svn
mailing list
|