Rivet analyses referenceOPAL_2000_I474010Correlations for $\Lambda^0$ and $\bar\Lambda^0$ production at LEP1Experiment: OPAL (LEP) Inspire ID: 474010 Status: VALIDATED Authors:
Beam energies: ANY Run details:
Measurement of the correlations between the production of $\Lambda^0$ and $\bar\Lambda^0$ baryons at LEP1, which is useful for testing models of baryon production. Source code: OPAL_2000_I474010.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/Thrust.hh"
5#include "Rivet/Projections/FastJets.hh"
6#include "Rivet/Projections/UnstableParticles.hh"
7
8
9namespace {
10
11 unsigned int factorial(unsigned int n) {
12 if(n<=1) return 0;
13 unsigned int out=1;
14 for (unsigned int i=n;i>0;--i) out *= i;
15 return out;
16 }
17
18 double rapidityAxis(Rivet::Vector4 momentum, Rivet::Vector3 axis) {
19 // assume axis is unit vector
20 double plong = momentum.vector3().dot(axis);
21 return 0.5*log((momentum.t()+plong)/(momentum.t()-plong));
22 }
23
24 double cosThetaStar(const Rivet::Particle & p1, const Rivet::Particle & p2, Rivet::Vector3 axis) {
25 Rivet::FourMomentum psum = p1.momentum()+p2.momentum();
26 Rivet::LorentzTransform trans(Rivet::LorentzTransform::mkObjTransformFromBeta(-psum.betaVec()));
27 Rivet::FourMomentum m1 = trans.transform(p1.momentum());
28 Rivet::FourVector a4;
29 a4.setX(axis.x());
30 a4.setY(axis.y());
31 a4.setZ(axis.z());
32 a4.setT(1.);
33 a4 = trans.transform(a4);
34 return fabs(m1.vector3().dot(a4.vector3())/m1.vector3().mod()/a4.vector3().mod());
35 }
36
37}
38
39namespace Rivet {
40
41
42 /// @brief lambda anti-lambda correlations
43 class OPAL_2000_I474010 : public Analysis {
44 public:
45
46 /// Constructor
47 RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_2000_I474010);
48
49
50 /// @name Analysis methods
51 /// @{
52
53 /// Book histograms and initialise projections before the run
54 void init() {
55 // projections
56 const FinalState fs;
57 declare(fs, "FS");
58 FastJets durhamjets(fs, JetAlg::DURHAM, 0.7);
59 durhamjets.useInvisibles(JetInvisibles::ALL);
60 declare(durhamjets, "DurhamJets");
61 const Thrust thrust(fs);
62 declare(thrust, "Thrust");
63 declare(UnstableParticles(), "UFS");
64 // histograms
65 book(_histLLBar , 1, 1, 1);
66 book(_histLL , 1, 1, 2);
67 book(_histLLBarCorr, 1, 1, 3);
68 book(_histLLBar_2jet , 2, 1, 1);
69 book(_histLL_2jet , 2, 1, 2);
70 book(_histLLBarCorr_2jet, 2, 1, 3);
71 book(_histLLBar_3jet , 3, 1, 1);
72 book(_histLL_3jet , 3, 1, 2);
73 book(_histLLBarCorr_3jet, 3, 1, 3);
74 book(_histdcosTheta, 4, 1, 1);
75 book(_histdy , 5, 1, 1);
76 // weights
77 book(_weight2Jet,"/TMP/W2Jet");
78 book(_weight3Jet,"/TMP/W3Jet");
79
80 }
81
82
83 /// Perform the per-event analysis
84 void analyze(const Event& event) {
85 // First, veto on leptonic events by requiring at least 4 charged FS particles
86 const FinalState& fs = apply<FinalState>(event, "FS");
87 const size_t numParticles = fs.particles().size();
88
89 // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
90 if (numParticles < 2) {
91 MSG_DEBUG("Failed leptonic event cut");
92 vetoEvent;
93 }
94 MSG_DEBUG("Passed leptonic event cut");
95
96 const Thrust& thrust = apply<Thrust>(event, "Thrust");
97 Vector3 axis = thrust.thrustAxis();
98
99 // Final state of unstable particles to get particle spectra
100 const UnstableParticles & ufs = apply<UnstableParticles>(event, "UFS");
101
102 Particles lambda,lambdaBar;
103 for (const Particle& p : ufs.particles()) {
104 const int id = p.pid();
105 if(id==3122) lambda.push_back(p);
106 else if(id==-3122) lambdaBar.push_back(p);
107 }
108
109 // get the number of jets
110 const FastJets& durjet = apply<FastJets>(event, "DurhamJets");
111
112 unsigned int njet = durjet.clusterSeq()->n_exclusive_jets_ycut(0.005);
113 if(njet==2) _weight2Jet->fill();
114 else if(njet==3) _weight3Jet->fill();
115
116 // need a pair of baryons
117 if(lambda.size()+lambdaBar.size()<2) vetoEvent;
118
119 // multiplicities
120 int nsame = (factorial(lambda.size())+factorial(lambdaBar.size()))/2;
121 int ndiff = lambda.size()*lambdaBar.size();
122 _histLLBar ->fill(Ecm, ndiff);
123 _histLL ->fill(Ecm, nsame);
124 _histLLBarCorr->fill(Ecm, (ndiff-nsame));
125
126 // multiplicities for different numbers of jets
127 if(njet==2) {
128 _histLLBar_2jet ->fill(Ecm, ndiff);
129 _histLL_2jet ->fill(Ecm, nsame);
130 _histLLBarCorr_2jet->fill(Ecm, (ndiff-nsame));
131 }
132 else if(njet==3) {
133 _histLLBar_3jet ->fill(Ecm, ndiff);
134 _histLL_3jet ->fill(Ecm, nsame);
135 _histLLBarCorr_3jet->fill(Ecm, (ndiff-nsame));
136 }
137 // uncorrelated lambda
138 for (unsigned int ix=0;ix<lambda.size();++ix) {
139 double y1 = rapidityAxis(lambda[ix].momentum(),axis);
140 for (unsigned int iy=ix+1;iy<lambda.size();++iy) {
141 double y2 = rapidityAxis(lambda[iy].momentum(),axis);
142 if(njet==2) _histdy->fill(fabs(y1-y2),-1.0);
143 double ctheta = cosThetaStar(lambda[ix],lambda[iy],axis);
144 _histdcosTheta->fill(ctheta,-1.0);
145 }
146 }
147 // uncorrelated lambdabar
148 for (unsigned int ix=0;ix<lambdaBar.size();++ix) {
149 double y1 = rapidityAxis(lambdaBar[ix].momentum(),axis);
150 for(unsigned int iy=ix+1;iy<lambdaBar.size();++iy) {
151 double y2 = rapidityAxis(lambdaBar[iy].momentum(),axis);
152 if(njet==2) _histdy->fill(fabs(y1-y2),-1.0);
153 double ctheta = cosThetaStar(lambdaBar[ix],lambdaBar[iy],axis);
154 _histdcosTheta->fill(ctheta,-1.0);
155 }
156 }
157 // all lambda lambdabar
158 for (unsigned int ix=0;ix<lambda.size();++ix) {
159 double y1 = rapidityAxis(lambda[ix].momentum(),axis);
160 for(unsigned int iy=0;iy<lambdaBar.size();++iy) {
161 double y2 = rapidityAxis(lambdaBar[iy].momentum(),axis);
162 if(njet==2) _histdy->fill(fabs(y1-y2));
163 double ctheta = cosThetaStar(lambda[ix],lambdaBar[iy],axis);
164 _histdcosTheta->fill(ctheta);
165 }
166 }
167 }
168
169
170 /// Normalise histograms etc., after the run
171 void finalize() {
172 scale(_histLLBar,1./sumOfWeights());
173 scale(_histLL,1./sumOfWeights());
174 scale(_histLLBarCorr,1./sumOfWeights());
175 scale(_histLLBar_2jet,1./ *_weight2Jet);
176 scale(_histLL_2jet,1./ *_weight2Jet);
177 scale(_histLLBarCorr_2jet,1./ *_weight2Jet);
178 scale(_histLLBar_3jet,1./ *_weight3Jet);
179 scale(_histLL_3jet,1./ *_weight3Jet);
180 scale(_histLLBarCorr_3jet,1./ *_weight3Jet);
181 normalize(_histdcosTheta, 1.);
182 normalize(_histdy , 1.);
183 }
184
185 /// @}
186
187
188 /// @name Histograms
189 /// @{
190
191 CounterPtr _weight2Jet;
192 CounterPtr _weight3Jet;
193
194 BinnedHistoPtr<string> _histLLBar, _histLL, _histLLBarCorr;
195 BinnedHistoPtr<string> _histLLBar_2jet, _histLL_2jet, _histLLBarCorr_2jet;
196 BinnedHistoPtr<string> _histLLBar_3jet, _histLL_3jet, _histLLBarCorr_3jet;
197
198 Histo1DPtr _histdcosTheta;
199 Histo1DPtr _histdy;
200
201 const string Ecm = "91.2";
202 /// @}
203
204
205 };
206
207
208 RIVET_DECLARE_PLUGIN(OPAL_2000_I474010);
209
210
211}
|