Rivet analyses referenceBABAR_2003_I593379Measurement of inclusive charmonium productionExperiment: BaBar (PEP-II) Inspire ID: 593379 Status: VALIDATED Authors:
Beam energies: (3.5, 8.0) GeV Run details:
Measurement of J/ψ, ψ′, χc1 and χc2 production using a data sample corresponding to an integrated luminosity of 20.3 fb−1 collected with the BABAR detector at the SLAC PEP-II electron-positron storage ring operating at a centre-of-mass energy near 10.58 GeV. Source code: BABAR_2003_I593379.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/UnstableParticles.hh"
4
5namespace Rivet {
6
7
8 /// @brief Babar charmonium spectra
9 /// @author Peter Richardson
10 class BABAR_2003_I593379 : public Analysis {
11 public:
12
13 BABAR_2003_I593379()
14 : Analysis("BABAR_2003_I593379")
15 { }
16
17
18 void analyze(const Event& e) {
19 // Find the charmonia
20 Particles upsilons;
21 // First in unstable final state
22 const UnstableParticles& ufs = apply<UnstableParticles>(e, "UFS");
23 for (const Particle& p : ufs.particles())
24 if (p.pid() == 300553) upsilons.push_back(p);
25 // Then in whole event if fails
26 if (upsilons.empty()) {
27 for(ConstGenParticlePtr p: HepMCUtils::particles(e.genEvent())) {
28 if (p->pdg_id() != 300553) continue;
29 ConstGenVertexPtr pv = p->production_vertex();
30 bool passed = true;
31 if (pv) {
32 for(ConstGenParticlePtr pp: HepMCUtils::particles(pv, Relatives::PARENTS)){
33 if ( p->pdg_id() == pp->pdg_id() ) {
34 passed = false;
35 break;
36 }
37 }
38 }
39 if (passed) upsilons.push_back(Particle(p));
40 }
41 }
42
43 // Find upsilons
44 for (const Particle& p : upsilons) {
45 _weightSum->fill();
46 // Find the charmonium resonances
47 /// @todo Use Rivet::Particles
48 vector<ConstGenParticlePtr> allJpsi, primaryJpsi, Psiprime, all_chi_c1, all_chi_c2, primary_chi_c1, primary_chi_c2;
49 findDecayProducts(p.genParticle(), allJpsi, primaryJpsi, Psiprime,
50 all_chi_c1, all_chi_c2, primary_chi_c1, primary_chi_c2);
51 const LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(p.mom().betaVec());
52 for (size_t i = 0; i < allJpsi.size(); i++) {
53 const double pcm = cms_boost.transform(FourMomentum(allJpsi[i]->momentum())).p();
54 _hist_all_Jpsi->fill(pcm);
55 }
56 _mult_JPsi->fill(double(allJpsi.size()));
57 for (size_t i = 0; i < primaryJpsi.size(); i++) {
58 const double pcm = cms_boost.transform(FourMomentum(primaryJpsi[i]->momentum())).p();
59 _hist_primary_Jpsi->fill(pcm);
60 }
61 _mult_JPsi_direct->fill(double(primaryJpsi.size()));
62 for (size_t i=0; i<Psiprime.size(); i++) {
63 const double pcm = cms_boost.transform(FourMomentum(Psiprime[i]->momentum())).p();
64 _hist_Psi_prime->fill(pcm);
65 }
66 _mult_Psi2S->fill(double(Psiprime.size()));
67 for (size_t i = 0; i < all_chi_c1.size(); i++) {
68 const double pcm = cms_boost.transform(FourMomentum(all_chi_c1[i]->momentum())).p();
69 _hist_chi_c1->fill(pcm);
70 }
71 _mult_chi_c1->fill(double(all_chi_c1.size()));
72 _mult_chi_c1_direct->fill(double(primary_chi_c1.size()));
73 for (size_t i = 0; i < all_chi_c2.size(); i++) {
74 const double pcm = cms_boost.transform(FourMomentum(all_chi_c2[i]->momentum())).p();
75 _hist_chi_c2->fill(pcm);
76 }
77 _mult_chi_c2->fill(double(all_chi_c2.size()));
78 _mult_chi_c2_direct->fill(double(primary_chi_c2.size()));
79 }
80 } // analyze
81
82
83 void finalize() {
84 scale(_hist_all_Jpsi , 0.5*0.1 / *_weightSum);
85 scale(_hist_chi_c1 , 0.5*0.1 / *_weightSum);
86 scale(_hist_chi_c2 , 0.5*0.1 / *_weightSum);
87 scale(_hist_Psi_prime , 0.5*0.1 / *_weightSum);
88 scale(_hist_primary_Jpsi , 0.5*0.1 / *_weightSum);
89 scale(_mult_JPsi , 0.5*100. / *_weightSum);
90 scale(_mult_JPsi_direct , 0.5*100. / *_weightSum);
91 scale(_mult_chi_c1 , 0.5*100. / *_weightSum);
92 scale(_mult_chi_c1_direct, 0.5*100. / *_weightSum);
93 scale(_mult_chi_c2 , 0.5*100. / *_weightSum);
94 scale(_mult_chi_c2_direct, 0.5*100. / *_weightSum);
95 scale(_mult_Psi2S , 0.5*100. / *_weightSum);
96 } // finalize
97
98
99 void init() {
100 declare(UnstableParticles(), "UFS");
101
102 book(_mult_JPsi ,1, 1, 1);
103 book(_mult_JPsi_direct ,1, 1, 2);
104 book(_mult_chi_c1 ,1, 1, 3);
105 book(_mult_chi_c1_direct ,1, 1, 4);
106 book(_mult_chi_c2 ,1, 1, 5);
107 book(_mult_chi_c2_direct ,1, 1, 6);
108 book(_mult_Psi2S ,1, 1, 7);
109 book(_hist_all_Jpsi ,2, 1, 1);
110 book(_hist_chi_c1 ,3, 1, 1);
111 book(_hist_chi_c2 ,3, 1, 2);
112 book(_hist_Psi_prime ,4, 1, 1);
113 book(_hist_primary_Jpsi ,6, 1, 1);
114
115 book(_weightSum, "TMP/weightSum");
116 } // init
117
118 private:
119
120 /// @{
121 // count of weights
122 CounterPtr _weightSum;
123 /// Histograms
124 Histo1DPtr _hist_all_Jpsi;
125 Histo1DPtr _hist_chi_c1;
126 Histo1DPtr _hist_chi_c2;
127 Histo1DPtr _hist_Psi_prime;
128 Histo1DPtr _hist_primary_Jpsi;
129
130 CounterPtr _mult_JPsi;
131 CounterPtr _mult_JPsi_direct;
132 CounterPtr _mult_chi_c1;
133 CounterPtr _mult_chi_c1_direct;
134 CounterPtr _mult_chi_c2;
135 CounterPtr _mult_chi_c2_direct;
136 CounterPtr _mult_Psi2S;
137 /// @}
138
139 void findDecayProducts(ConstGenParticlePtr p,
140 vector<ConstGenParticlePtr>& allJpsi,
141 vector<ConstGenParticlePtr>& primaryJpsi,
142 vector<ConstGenParticlePtr>& Psiprime,
143 vector<ConstGenParticlePtr>& all_chi_c1, vector<ConstGenParticlePtr>& all_chi_c2,
144 vector<ConstGenParticlePtr>& primary_chi_c1, vector<ConstGenParticlePtr>& primary_chi_c2) {
145 ConstGenVertexPtr dv = p->end_vertex();
146 bool isOnium = false;
147 /// @todo Use better looping
148 for (ConstGenParticlePtr pp: HepMCUtils::particles(dv, Relatives::PARENTS)){
149 int id = pp->pdg_id();
150 id = id%1000;
151 id -= id%10;
152 id /= 10;
153 if (id==44) isOnium = true;
154 }
155 /// @todo Use better looping
156 for (ConstGenParticlePtr pp: HepMCUtils::particles(dv, Relatives::CHILDREN)){
157 int id = pp->pdg_id();
158 if (id==100443) {
159 Psiprime.push_back(pp);
160 }
161 else if (id==20443) {
162 all_chi_c1.push_back(pp);
163 if (!isOnium) primary_chi_c1.push_back(pp);
164 }
165 else if (id==445) {
166 all_chi_c2.push_back(pp);
167 if (!isOnium) primary_chi_c2.push_back(pp);
168 }
169 else if (id==443) {
170 allJpsi.push_back(pp);
171 if (!isOnium) primaryJpsi.push_back(pp);
172 }
173 if (pp->end_vertex()) {
174 findDecayProducts(pp, allJpsi, primaryJpsi, Psiprime, all_chi_c1, all_chi_c2, primary_chi_c1, primary_chi_c2);
175 }
176 }
177 }
178
179 };
180
181
182 RIVET_DECLARE_PLUGIN(BABAR_2003_I593379);
183
184}
|
