Rivet analyses referenceATLAS_2011_I944826$K_S^0$ and $\Lambda$ production at 0.9 and 7 TeV with ATLASExperiment: ATLAS (LHC) Inspire ID: 944826 Status: VALIDATED Authors:
Beam energies: (450.0, 450.0); (3500.0, 3500.0) GeV Run details:
The production of $K_S$ and $\Lambda$ hadrons is studied in inelastic $pp$ collisions at $\sqrt{s} = 0.9$ and 7 TeV collected with the ATLAS detector at the LHC using a minimum-bias trigger. The observed distributions of transverse momentum, rapidity, and multiplicity are corrected to hadron level in a model-independent way within well defined phase-space regions. The distribution of the production ratio of $\bar{\Lambda}$ to $\Lambda$ baryons is also measured. The results are compared with various Monte Carlo simulation models. Although most of these models agree with data to within 15 percent in the $K_S$ distributions, substantial disagreements with data are found in the $\Lambda$ distributions of transverse momentum. Source code: ATLAS_2011_I944826.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/ChargedFinalState.hh"
4#include "Rivet/Projections/IdentifiedFinalState.hh"
5#include "Rivet/Projections/UnstableParticles.hh"
6
7namespace Rivet {
8
9
10 class ATLAS_2011_I944826 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_I944826);
15
16 /// Book histograms and initialise projections before the run
17 void init() {
18
19 UnstableParticles ufs(Cuts::pT > 100*MeV);
20 declare(ufs, "UFS");
21
22 ChargedFinalState mbts(Cuts::absetaIn(2.09, 3.84));
23 declare(mbts, "MBTS");
24
25 IdentifiedFinalState nstable(Cuts::abseta < 2.5 && Cuts::pT >= 100*MeV);
26 nstable.acceptIdPair(PID::ELECTRON)
27 .acceptIdPair(PID::MUON)
28 .acceptIdPair(PID::PIPLUS)
29 .acceptIdPair(PID::KPLUS)
30 .acceptIdPair(PID::PROTON);
31 declare(nstable, "nstable");
32
33
34 for (double eVal : allowedEnergies()) {
35
36 const string en = toString(int(eVal));
37 if (isCompatibleWithSqrtS(eVal)) _sqs = en;
38
39 bool is900(en == "900");
40 size_t offset = is900? 3 : 0;
41 book(_h[en+"Ks_pT"], 1+offset, 1, 1);
42 book(_h[en+"Ks_y"], 2+offset, 1, 1);
43 book(_h[en+"Ks_mult"], 3+offset, 1, 1);
44 book(_h[en+"L_pT"], 7+offset, 1, 1);
45 book(_h[en+"L_y"], 8+offset, 1, 1);
46 book(_h[en+"L_mult"], 9+offset, 1, 1);
47 if (is900) offset = 2;
48 book(_e[en+"v_y"], 13+offset, 1, 1);
49 book(_e[en+"v_pT"], 14+offset, 1, 1);
50 //
51 book(_h[en+"lambda_v_y"], "TMP/lambda_v_y_"+en, is900? 5 : 10, 0.0, 2.5);
52 book(_h[en+"lambdabar_v_y"], "TMP/lambdabar_v_y_"+en, is900? 5 : 10, 0.0, 2.5);
53 book(_h[en+"lambda_v_pT"], "TMP/lambda_v_pT_"+en, is900? 8 : 18, 0.5, is900? 3.7 : 4.1);
54 book(_h[en+"lambdabar_v_pT"], "TMP/lambdabar_v_pT_"+en, is900? 8 : 18, 0.5, is900? 3.7 : 4.1);
55 }
56 if (_sqs == "" && !merging()) {
57 throw BeamError("Invalid beam energy for " + name() + "\n");
58 }
59 }
60
61
62 // This function is required to impose the flight time cuts on Kaons and Lambdas
63 double getPerpFlightDistance(const Rivet::Particle& p) {
64 ConstGenParticlePtr genp = p.genParticle();
65 ConstGenVertexPtr prodV = genp->production_vertex();
66 ConstGenVertexPtr decV = genp->end_vertex();
67 RivetHepMC::FourVector prodPos = prodV->position();
68 if (decV) {
69 const RivetHepMC::FourVector decPos = decV->position();
70 double dy = prodPos.y() - decPos.y();
71 double dx = prodPos.x() - decPos.x();
72 return add_quad(dx, dy);
73 }
74 return numeric_limits<double>::max();
75 }
76
77
78 bool daughtersSurviveCuts(const Rivet::Particle& p) {
79 // We require the Kshort or Lambda to decay into two charged
80 // particles with at least pT = 100 MeV inside acceptance region
81 ConstGenParticlePtr genp = p.genParticle();
82 ConstGenVertexPtr decV = genp->end_vertex();
83 bool decision = true;
84
85 if (!decV) return false;
86 if (HepMCUtils::particles(decV, Relatives::CHILDREN).size() == 2) {
87 std::vector<double> pTs;
88 std::vector<int> charges;
89 std::vector<double> etas;
90 for(ConstGenParticlePtr gp: HepMCUtils::particles(decV, Relatives::CHILDREN)) {
91 pTs.push_back(gp->momentum().perp());
92 etas.push_back(fabs(gp->momentum().eta()));
93 charges.push_back( Rivet::PID::charge3(gp->pdg_id()) );
94 // gp->print();
95 }
96 if ( (pTs[0]/Rivet::GeV < 0.1) || (pTs[1]/Rivet::GeV < 0.1) ) {
97 decision = false;
98 MSG_DEBUG("Failed pT cut: " << pTs[0]/Rivet::GeV << " " << pTs[1]/Rivet::GeV);
99 }
100 if ( etas[0] > 2.5 || etas[1] > 2.5 ) {
101 decision = false;
102 MSG_DEBUG("Failed eta cut: " << etas[0] << " " << etas[1]);
103 }
104 if ( charges[0] * charges[1] >= 0 ) {
105 decision = false;
106 MSG_DEBUG("Failed opposite charge cut: " << charges[0] << " " << charges[1]);
107 }
108 }
109 else {
110 decision = false;
111 MSG_DEBUG("Failed nDaughters cut: " << HepMCUtils::particles(decV, Relatives::CHILDREN).size());
112 }
113
114 return decision;
115 }
116
117 /// Perform the per-event analysis
118 void analyze(const Event& event) {
119
120 // ATLAS MBTS trigger requirement of at least one hit in either hemisphere
121 if (apply<FinalState>(event, "MBTS").size() < 1) {
122 MSG_DEBUG("Failed trigger cut");
123 vetoEvent;
124 }
125
126 // Veto event also when we find less than 2 particles in the acceptance region of type 211,2212,11,13,321
127 if (apply<FinalState>(event, "nstable").size() < 2) {
128 MSG_DEBUG("Failed stable particle cut");
129 vetoEvent;
130 }
131
132 // This ufs holds all the Kaons and Lambdas
133 const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
134
135 // Some conters
136 int n_KS0 = 0, n_LAMBDA = 0;
137
138 // Particle loop
139 for (const Particle& p : ufs.particles()) {
140
141 // General particle quantities
142 const double pT = p.pT();
143 const double y = p.rapidity();
144 const PdgId apid = p.abspid();
145
146 double flightd = 0.0;
147
148 // Look for Kaons, Lambdas
149 switch (apid) {
150
151 case PID::K0S:
152 flightd = getPerpFlightDistance(p);
153 if (!inRange(flightd/mm, 4., 450.) ) {
154 MSG_DEBUG("Kaon failed flight distance cut:" << flightd);
155 break;
156 }
157 if (daughtersSurviveCuts(p) ) {
158 _h[_sqs+"Ks_y"] ->fill(y);
159 _h[_sqs+"Ks_pT"]->fill(pT/GeV);
160 ++n_KS0;
161 }
162 break;
163
164 case PID::LAMBDA:
165 if (pT < 0.5*GeV) { // Lambdas have an additional pT cut of 500 MeV
166 MSG_DEBUG("Lambda failed pT cut:" << pT/GeV << " GeV");
167 break;
168 }
169 flightd = getPerpFlightDistance(p);
170 if (!inRange(flightd/mm, 17., 450.)) {
171 MSG_DEBUG("Lambda failed flight distance cut:" << flightd/mm << " mm");
172 break;
173 }
174 if ( daughtersSurviveCuts(p) ) {
175 if (p.pid() == PID::LAMBDA) {
176 _h[_sqs+"lambda_v_y"]->fill(fabs(y));
177 _h[_sqs+"lambda_v_pT"]->fill(pT/GeV);
178 _h[_sqs+"L_y"]->fill(y);
179 _h[_sqs+"L_pT"]->fill(pT/GeV);
180 ++n_LAMBDA;
181 }
182 else if (p.pid() == -PID::LAMBDA) {
183 _h[_sqs+"lambdabar_v_y"]->fill(fabs(y));
184 _h[_sqs+"lambdabar_v_pT"]->fill(pT/GeV);
185 }
186 }
187 break;
188 }
189 }
190
191 // Fill multiplicity histos
192 _h[_sqs+"Ks_mult"]->fill(n_KS0);
193 _h[_sqs+"L_mult"]->fill(n_LAMBDA);
194 }
195
196
197
198 /// Normalise histograms etc., after the run
199 void finalize() {
200 // Division of histograms to obtain lambda_bar/lambda ratios
201 for (double eVal : allowedEnergies()) {
202 const string en = toString(int(eVal));
203
204 if (_h[en+"L_pT"]->sumW()) scale(_h[en+"L_y"], 1.0 / _h[en+"L_pT"]->sumW());
205 divide(_h[en+"lambdabar_v_y"], _h[en+"lambda_v_y"], _e[en+"v_y"]);
206 divide(_h[en+"lambdabar_v_pT"], _h[en+"lambda_v_pT"], _e[en+"v_pT"]);
207 }
208 normalize(_h);
209
210 }
211
212
213 private:
214
215 /// @name Persistent histograms
216 /// @{
217 map<string,Histo1DPtr> _h;
218 map<string,Estimate1DPtr> _e;
219
220 string _sqs = "";
221 /// @}
222
223 };
224
225
226 RIVET_DECLARE_PLUGIN(ATLAS_2011_I944826);
227
228}
|