1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
| // -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/Beam.hh"
#include "Rivet/Projections/UnstableParticles.hh"
namespace Rivet {
/// @brief BABAR Xi_c baryons from fragmentation
///
/// @author Peter Richardson
class BABAR_2005_S6181155 : public Analysis {
public:
RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2005_S6181155);
void init() {
declare(Beam(), "Beams");
declare(UnstableParticles(), "UFS");
book(_histOnResonanceA, 1,1,1);
book(_histOnResonanceB, 2,1,1);
book(_histOffResonance, 2,1,2);
book(_sigma, 3,1,1);
}
void analyze(const Event& e) {
// Loop through unstable FS particles and look for charmed mesons/baryons
const UnstableParticles& ufs = apply<UnstableParticles>(e, "UFS");
const Beam beamproj = apply<Beam>(e, "Beams");
const ParticlePair& beams = beamproj.beams();
const FourMomentum mom_tot = beams.first.momentum() + beams.second.momentum();
const LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(mom_tot.betaVec());
const double s = sqr(beamproj.sqrtS());
const bool onresonance = fuzzyEquals(beamproj.sqrtS()/GeV, 10.58, 2E-3);
for (const Particle& p : ufs.particles()) {
// 3-momentum in CMS frame
const double mom = cms_boost.transform(p.momentum()).vector3().mod();
// Only looking at Xi_c^0
if (p.abspid() != 4132 ) continue;
MSG_DEBUG("mom = " << mom);
// off-resonance cross section
if (checkDecay(p.genParticle())) {
if (onresonance) {
_histOnResonanceA->fill(mom);
_histOnResonanceB->fill(mom);
}
else {
_histOffResonance->fill(mom,s/sqr(10.58));
_sigma->fill(10.58);
}
}
}
}
void finalize() {
scale(_histOnResonanceA, crossSection()/femtobarn/sumOfWeights()*0.2);
scale(_histOnResonanceB, crossSection()/femtobarn/sumOfWeights()*0.45);
scale(_histOffResonance, crossSection()/femtobarn/sumOfWeights()*0.45);
scale(_sigma , crossSection()/femtobarn/sumOfWeights());
}
private:
/// @name Histograms
/// @{
Histo1DPtr _histOnResonanceA;
Histo1DPtr _histOnResonanceB;
Histo1DPtr _histOffResonance;
Histo1DPtr _sigma;
/// @}
bool checkDecay(ConstGenParticlePtr p) {
unsigned int nstable = 0, npip = 0, npim = 0;
unsigned int nXim = 0, nXip = 0;
findDecayProducts(p, nstable, npip, npim, nXip, nXim);
int id = p->pdg_id();
// Xi_c
if (id == 4132) {
if (nstable == 2 && nXim == 1 && npip == 1) return true;
}
else if (id == -4132) {
if (nstable == 2 && nXip == 1 && npim == 1) return true;
}
return false;
}
void findDecayProducts(ConstGenParticlePtr p,
unsigned int& nstable,
unsigned int& npip, unsigned int& npim,
unsigned int& nXip, unsigned int& nXim) {
ConstGenVertexPtr dv = p->end_vertex();
/// @todo Use better looping
for (ConstGenParticlePtr pp: HepMCUtils::particles(dv, Relatives::CHILDREN)){
int id = pp->pdg_id();
if (id==3312) {
++nXim;
++nstable;
} else if (id == -3312) {
++nXip;
++nstable;
} else if(id == 111 || id == 221) {
++nstable;
} else if (pp->end_vertex()) {
findDecayProducts(pp, nstable, npip, npim, nXip, nXim);
} else {
if (id != 22) ++nstable;
if (id == 211) ++npip;
else if(id == -211) ++npim;
}
}
}
};
RIVET_DECLARE_ALIASED_PLUGIN(BABAR_2005_S6181155, BABAR_2005_I679961);
}
|