Page MenuHomeHEPForge
Files F8310200

No OneTemporary
Actions

View Options

Index: trunk/src/Analyses/ARGUS_1993_S2669951.cc
===================================================================
--- trunk/src/Analyses/ARGUS_1993_S2669951.cc (revision 4288)
+++ trunk/src/Analyses/ARGUS_1993_S2669951.cc (revision 4289)
@@ -1,207 +1,205 @@
// -*- C++ -*-
#include <iostream>
#include "Rivet/Analysis.hh"
#include "Rivet/RivetYODA.hh"
#include "Rivet/Tools/ParticleIdUtils.hh"
#include "Rivet/Projections/Beam.hh"
#include "Rivet/Projections/UnstableFinalState.hh"
#include "Rivet/ParticleName.hh"
namespace Rivet {
/// @brief BELLE pi+/-, K+/- and proton/antiproton spectrum at Upsilon(4S)
/// @author Peter Richardson
class ARGUS_1993_S2669951 : public Analysis {
public:
ARGUS_1993_S2669951()
: Analysis("ARGUS_1993_S2669951"),
_count_etaPrime_highZ(2, 0.), _count_etaPrime_allZ(3, 0.), _count_f0(3, 0.),
_weightSum_cont(0.), _weightSum_Ups1(0.), _weightSum_Ups2(0.)
{ }
void init() {
addProjection(Beam(), "Beams");
addProjection(UnstableFinalState(), "UFS");
_hist_cont_f0 = bookHisto1D(2, 1, 1);
_hist_Ups1_f0 = bookHisto1D(3, 1, 1);
_hist_Ups2_f0 = bookHisto1D(4, 1, 1);
}
void analyze(const Event& e) {
const double weight = e.weight();
const Beam beamproj = applyProjection<Beam>(e, "Beams");
const double s = sqr(beamproj.sqrtS());
const double roots = sqrt(s);
const UnstableFinalState& ufs = applyProjection<UnstableFinalState>(e, "UFS");
// Find the upsilons
Particles upsilons;
// First in unstable final state
foreach (const Particle& p, ufs.particles())
if (p.pdgId() == 553 || p.pdgId() == 100553 ) upsilons.push_back(p);
// Then in whole event if fails
if (upsilons.empty()) {
foreach (GenParticle* p, Rivet::particles(e.genEvent())) {
if ( p->pdg_id() != 553 && p->pdg_id() != 100553 ) continue;
const GenVertex* pv = p->production_vertex();
bool passed = true;
if (pv) {
for (GenVertex::particles_in_const_iterator pp = pv->particles_in_const_begin(); pp != pv->particles_in_const_end() ; ++pp) {
if ( p->pdg_id() == (*pp)->pdg_id() ) {
passed = false;
break;
}
}
}
if (passed) upsilons.push_back(Particle(*p));
}
}
// Continuum
if (upsilons.empty()) {
_weightSum_cont += weight;
unsigned int nEtaA(0),nEtaB(0),nf0(0);
foreach (const Particle& p, ufs.particles()) {
int id = abs(p.pdgId());
double xp = 2.*p.momentum().t()/roots;
double beta = p.momentum().vector3().mod() / p.momentum().t();
if (id == 9010221) {
_hist_cont_f0->fill(xp, weight/beta);
++nf0;
} else if (id == 331) {
if (xp > 0.35) ++nEtaA;
++nEtaB;
}
}
_count_f0[2] += nf0*weight;
_count_etaPrime_highZ[1] += nEtaA*weight;
_count_etaPrime_allZ[2] += nEtaB*weight;
}
else {
// find an upsilons
foreach (const Particle& ups, upsilons) {
int parentId = ups.pdgId();
((parentId == 553) ? _weightSum_Ups1 : _weightSum_Ups2) += weight;
Particles unstable;
// find the decay products we want
findDecayProducts(ups.genParticle(), unstable);
LorentzTransform cms_boost;
if (ups.momentum().vector3().mod() > 1*MeV)
cms_boost = LorentzTransform(-ups.momentum().boostVector());
double mass = ups.momentum().mass();
unsigned int nEtaA(0), nEtaB(0), nf0(0);
foreach(const Particle& p , unstable) {
const int id = abs(p.pdgId());
FourMomentum p2 = cms_boost.transform(p.momentum());
double xp = 2.*p2.t()/mass;
double beta = p2.vector3().mod()/p2.t();
if (id == 9010221) {
((parentId == 553) ? _hist_Ups1_f0 : _hist_Ups2_f0)->fill(xp, weight/beta);
++nf0;
} else if (id == 331) {
if (xp > 0.35) ++nEtaA;
++nEtaB;
}
}
if (parentId == 553) {
_count_f0[0] += nf0*weight;
_count_etaPrime_highZ[0] += nEtaA*weight;
_count_etaPrime_allZ[0] += nEtaB*weight;
} else {
_count_f0[1] += nf0*weight;
_count_etaPrime_allZ[1] += nEtaB*weight;
}
}
}
}
void finalize() {
// @todo YODA
- Scatter2DPtr mult_etaPrime_highZ = bookScatter2D(1, 1, 1);
- for (size_t i = 0; i < mult_etaPrime_highZ->numPoints(); ++i) {
- Point2D& p = mult_etaPrime_highZ->point(i);
+ foreach (Point2D& p, bookScatter2D(1, 1, 1)->points()) {
if (fuzzyEquals(9.905, p.x(), 1e-3) && _weightSum_cont > 0) {
p.setY(_count_etaPrime_highZ[1] / _weightSum_cont);
} else if (fuzzyEquals( 9.46, p.x(), 1e-3) && _weightSum_Ups1 > 0) {
p.setY(_count_etaPrime_highZ[0] / _weightSum_Ups1);
}
}
//AIDA::IDataPointSet * mult_etaPrime_allZ = bookDataPointSet( 1,1,2);
//for (int i = 0; i < mult_etaPrime_allZ->size(); ++i) {
// if ( fuzzyEquals( 9.905, mult_etaPrime_allZ->point(i)->coordinate(0)->value(), 0.001) &&
// _weightSum_cont>0.) {
// mult_etaPrime_allZ->point(i)->coordinate(1)->setValue( _count_etaPrime_allZ[2]/_weightSum_cont);
// }
// else if ( fuzzyEquals( 9.46, mult_etaPrime_allZ->point(i)->coordinate(0)->value(), 0.001) &&
// _weightSum_Ups1>0.) {
// mult_etaPrime_allZ->point(i)->coordinate(1)->setValue(_count_etaPrime_allZ[0]/_weightSum_Ups1);
// }
// else if ( fuzzyEquals( 10.02, mult_etaPrime_allZ->point(i)->coordinate(0)->value(), 0.001) &&
// _weightSum_Ups2>0.) {
// mult_etaPrime_allZ->point(i)->coordinate(1)->setValue(_count_etaPrime_allZ[1]/_weightSum_Ups2);
// }
//}
//AIDA::IDataPointSet * mult_f0 = bookDataPointSet( 5,1,1);
//for (int i = 0; i < mult_f0->size(); ++i) {
// if ( fuzzyEquals( 10.45, mult_f0->point(i)->coordinate(0)->value(), 0.001) &&
// _weightSum_cont>0.) {
// mult_f0->point(i)->coordinate(1)->setValue( _count_f0[2]/_weightSum_cont);
// }
// else if ( fuzzyEquals( 9.46, mult_f0->point(i)->coordinate(0)->value(), 0.001) &&
// _weightSum_Ups1>0.) {
// mult_f0->point(i)->coordinate(1)->setValue(_count_f0[0]/_weightSum_Ups1);
// }
// else if ( fuzzyEquals( 10.02, mult_f0->point(i)->coordinate(0)->value(), 0.001) &&
// _weightSum_Ups2>0.) {
// mult_f0->point(i)->coordinate(1)->setValue(_count_f0[1]/_weightSum_Ups2);
// }
//}
if (_weightSum_cont > 0.) scale(_hist_cont_f0, 1./_weightSum_cont);
if (_weightSum_Ups1 > 0.) scale(_hist_Ups1_f0, 1./_weightSum_Ups1);
if (_weightSum_Ups2 > 0.) scale(_hist_Ups2_f0, 1./_weightSum_Ups2);
}
private:
/// @name Counters
//@{
vector<double> _count_etaPrime_highZ, _count_etaPrime_allZ, _count_f0;
double _weightSum_cont,_weightSum_Ups1,_weightSum_Ups2;
//@}
/// Histos
Histo1DPtr _hist_cont_f0, _hist_Ups1_f0, _hist_Ups2_f0;
/// Recursively walk the HepMC tree to find decay products of @a p
void findDecayProducts(const GenParticle* p, Particles& unstable) {
const GenVertex* dv = p->end_vertex();
/// @todo Use better looping
for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin(); pp != dv->particles_out_const_end(); ++pp) {
const int id = abs((*pp)->pdg_id());
if (id == 331 || id == 9010221) unstable.push_back(Particle(*pp));
else if ((*pp)->end_vertex()) findDecayProducts(*pp, unstable);
}
}
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(ARGUS_1993_S2669951);
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 5:40 PM (12 h, 12 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023676
Default Alt Text
(8 KB)

Event Timeline