Rivet analyses referenceATLAS_2011_I944826K0S and Λ 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 KS and Λ hadrons is studied in inelastic pp collisions at √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 ˉΛ to Λ 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 KS distributions, substantial disagreements with data are found in the Λ 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}
|
