Rivet analyses referenceATLAS_2017_I1624693Study of ordered hadron chains at 7 TeVExperiment: ATLAS (LHC) Inspire ID: 1624693 Status: VALIDATED Authors:
Beam energies: (3500.0, 3500.0) GeV Run details:
The analysis of the momentum difference between charged hadrons in high-energy proton-proton collisions is performed in order to study coherent particle production. The observed correlation pattern agrees with a model of a helical QCD string fragmenting into a chain of ground-state hadrons. A threshold momentum difference in the production of adjacent pairs of charged hadrons is observed, in agreement with model predictions. The presence of low-mass hadron chains also explains the emergence of charge-combination-dependent two-particle correlations commonly attributed to Bose-Einstein interference. The data sample consists of 190 $\mu\text{b}^{-1}$ of minimum-bias events collected with proton-proton collisions at a center-of-mass energy $\sqrt{s}$ = 7 TeV in the early low-luminosity data taking with the ATLAS detector at the LHC. Source code: ATLAS_2017_I1624693.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/ChargedFinalState.hh"
5
6namespace Rivet {
7
8 /// @brief Study of ordered hadron chains at 7 TeV
9 class ATLAS_2017_I1624693 : public Analysis {
10 public:
11
12 /// Constructor
13 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1624693);
14
15 /// @name Analysis methods
16 /// @{
17 struct usedX {
18
19 int locMin;
20 int locMax;
21 std::vector<std::pair<int,float> > chains;
22
23 // Constructor
24 usedX(int min, int max, int ic, float mass) {
25 locMin=min;
26 locMax=max;
27 chains.clear();
28 chains.push_back(std::pair<int,float>(ic,mass));
29 }
30
31 // Constructor
32 usedX(int min, int max) {
33 locMin=min;
34 locMax=max;
35 chains.clear();
36 }
37
38 void add(int jc, float mass) {
39
40 if (chains.size()) {
41 std::vector<std::pair<int,float> >::iterator it=chains.begin();
42 while ( it!=chains.end() && mass>(*it).second ) ++it;
43 chains.insert(it,std::pair<int,float>(jc,mass));
44 }
45 else {
46 chains.push_back(std::pair<int,float>(jc,mass));
47 }
48 }
49 };
50
51
52 /// Book histograms and initialise projections before the run
53 void init() {
54
55 /// @todo Initialise and register projections here
56 ChargedFinalState cfs((Cuts::etaIn(-2.5, 2.5) && Cuts::pT >= 0.1*GeV));
57 declare(cfs,"CFS");
58
59 // pion mass;
60 pim = 0.1396;
61
62 /// @todo Book histograms here, e.g.:
63 book(_DeltaQ , 1, 1, 1);
64 book(_Delta3h, 2, 1, 1);
65 book(_dalitz , 3, 1, 1);
66
67 // auxiliary
68 book(_h_nch, "_nch", 200, -0.5, 199.5);
69
70 }
71
72
73 /// Perform the per-event analysis
74 void analyze(const Event& event) {
75 //const double weight = event.weight();
76 bool match =false;
77
78 /// @todo Do the event by event analysis here
79 const ChargedFinalState& had = apply<ChargedFinalState>(event, "CFS");
80 Particles hs=had.particles();
81 int nch = hs.size();
82
83 if (nch < 3) return;
84
85 _h_nch->fill(1.*nch,1.);
86
87 for (unsigned int i=0; i < hs.size() - 1; ++i) {
88 for (unsigned int j=i+1; j < hs.size(); ++j) {
89 double q12 = qq(hs[i],hs[j],match);
90 if (match) _DeltaQ->fill(q12,-1.);
91 else _DeltaQ->fill(q12,1.);
92 }
93 }
94
95 // chain selection
96
97 std::vector<float> wchain;
98 std::vector< std::vector<unsigned int> > rchains;
99 std::vector< std::vector<float> > mchains;
100 wchain.clear(); rchains.clear(); mchains.clear();
101 for (unsigned int ip1 = 0; ip1< hs.size(); ++ip1 ) {
102 wchain.push_back(1.);
103 std::vector<unsigned int> cc(1,ip1);
104 std::vector<float> mc;
105
106 double qlmin=10000.; int ilmin=-1;
107 for (unsigned ip2 = 0; ip2 < hs.size(); ++ip2) {
108 if (ip2==ip1) continue;
109 double ql = qq(hs[ip1],hs[ip2],match);
110 if (!match) continue; // looking for closest like-sign match
111 if (ql <qlmin) { qlmin=ql; ilmin=ip2;}
112 }
113 if (ilmin<0) {
114 wchain.back()=0.;
115 mc.push_back(-1.);
116 }
117 else { // search for unlike-sign match
118 cc.push_back(ilmin);
119 mc.push_back(qlmin);
120 if (int(ip1)>ilmin && rchains[ilmin][1]==ip1) {
121 // std::cout <<"exclusive match:"<< std::endl;
122 wchain.back()=0.5; wchain[ilmin]=0.5;
123 }
124
125 double m3min=10000.; int ixmin=-1;
126 for (unsigned ip2 = 0; ip2< hs.size(); ++ip2) {
127 if (ip2==ip1 || int(ip2)==ilmin ) continue;
128 double qx = qq(hs[ip1],hs[ip2],match);
129 if (match) continue;
130 double qxl = qq(hs[ip2],hs[ilmin],match);
131 double m3 = sqrt(9*pim*pim+qxl*qxl+qlmin*qlmin+qx*qx);
132 if (m3 <m3min) { m3min=m3; ixmin=ip2;}
133 }
134
135 if (ixmin<0) {
136 wchain.back()=0.;
137 mc.push_back(-1.);
138 }
139 else {
140 cc.push_back(ixmin);
141 mc.push_back(m3min);
142 }
143 }
144 rchains.push_back(cc);
145 mchains.push_back(mc);
146 }
147
148 // cleanup: association rate for like-sign pairs should not exceed 2
149 std::vector<float> assoc(hs.size(),0.); // cache for association rate
150 std::vector<bool> accept(rchains.size(), false);
151 // loop over chains and accept lowest masses while watching the association rate
152 int inext = 0;
153 while ( inext>-1 ) {
154 inext = -1; float cMin = 100000.;
155 // find non-accepted chain with lowest Q_ls; dissolve chains if association count over 2
156 for (unsigned int ic=0; ic < rchains.size(); ++ic) {
157 if (rchains[ic].size() < 2) continue;
158 if (accept[ic]) continue;
159 if (mchains[ic][0] < cMin) { cMin = mchains[ic][0]; inext=ic; }
160 }
161 if (inext>-1 ) {
162 unsigned int cloc0 = rchains[inext][0];
163 unsigned int cloc1 = rchains[inext][1];
164 if ( (assoc[cloc0] + 1. <= 2.) && (assoc[cloc1] + 1. <= 2.) ) { // chain can be accepted
165 accept[inext]=true;
166 assoc[cloc0]+=1.;
167 assoc[cloc1]+=1.;
168 if (wchain[inext]==0.5) { // accept the identical chain, too
169 for (unsigned int ic=0; ic<hs.size(); ++ic) {
170 if (rchains[ic][0] == cloc1 && rchains[ic][1] == cloc0) {
171 accept[ic]=true;
172 break;
173 }
174 }
175 }
176 }
177 else if ( assoc[cloc0]>1 ) { // association count filled up, discard chain
178 accept[inext]=true;
179 wchain[inext]=0.;
180 }
181 else { // dissolve chain and find new association
182 unsigned int i1 = rchains[inext][0];
183 float mMn = 1000000.;
184 int ipn = -1;
185 for (unsigned int i2=0; i2<hs.size(); ++i2) {
186 if (i1 == i2) continue;
187 double m = qq(hs[i1],hs[i2],match);
188 if (!match) continue;
189 if (assoc[i2] > 1.) continue;
190 if (m > 0. && m <mMn ) { mMn = m; ipn = i2;}
191 }
192 if (ipn >= 0) {
193 rchains[inext][1]=ipn; mchains[inext][0]=mMn;
194 // resolve chain weight : by default, it is 1.
195 wchain[inext]=1.;
196 // check exclusivity of pairing
197 for (unsigned int ico=0; ico<hs.size(); ++ico) {
198 if (int(rchains[ico][0]) == ipn && rchains[ico][1] == i1) { // scale the contribution from both chains
199 wchain[ico]=0.5;
200 wchain[inext]=0.5;
201 }
202 }
203 // add 3.member
204 // continue with arbitrary match
205 int ipnn=-1; float mMnn = 10000.;
206 mMn = 1000000.;
207 for (unsigned int ij=0; ij < hs.size(); ++ij) {
208 rchains[inext].resize(2);
209 float q02 = qq(hs[i1],hs[ij],match);
210 if (match>0.) continue;
211 float q12 = qq(hs[ipn],hs[ij],match);
212 double m3 = sqrt(9*pim*pim+q02*q02+mMn*mMn+q12*q12);
213 if (m3>0. && m3 <mMnn ) { mMnn = m3; ipnn = ij; }
214 }
215 if (ipnn>=0) { rchains[inext].push_back(ipnn); rchains[inext][2]=ipnn; mchains[inext][1]=mMnn; }
216 else {accept[inext]=true; wchain[inext]=0.;}
217 }
218 else { // chain not recovered
219 wchain[inext]=0.;
220 accept[inext]=true;
221 }
222 }
223 }
224 } // end loop over chains
225
226 // cleanup: association rate for unlike-sign pairs
227 // third member verification
228 std::vector<bool> accept3(rchains.size(),false);
229 // watch unlike-sign combinations used
230 std::vector<usedX> used;
231 // loop over chains and accept lowest masses while watching the association rate
232 inext = 0;
233 while ( inext>-1 ) {
234 inext = -1; float cMin = 100000.;
235 // find non-accepted chain with lowest mass; dissolve chains if association count over 3
236 for (unsigned int ic=0; ic < rchains.size(); ++ic) {
237 if (rchains[ic].size() < 3 || !wchain[ic] || !accept[ic]) continue;
238 if (accept3[ic]) continue;
239 if (mchains[ic][1]<cMin) { cMin = mchains[ic][1]; inext=ic; }
240 }
241 // check association counts
242 if (inext>-1 ) {
243 unsigned int cloc0 = rchains[inext][0];
244 unsigned int cloc1 = rchains[inext][1];
245 unsigned int cloc2 = rchains[inext][2];
246
247 // map use of unlike sign pairs
248 int iu0 = -1; float w0=0.;
249 for (unsigned int iu=0; iu<used.size(); ++iu) {
250 if (fmin(cloc0,cloc2)==used[iu].locMin && fmax(cloc0,cloc2)==used[iu].locMax ) {
251 iu0=iu;
252 if (used[iu].chains.size() > 0)
253 for (unsigned int iw=0; iw<used[iu].chains.size(); ++iw) w0+=wchain[used[iu].chains[iw].first];
254 //used[iu].add(i1,mch[1]);
255 break;
256 }
257 }
258 if (iu0<0) { used.push_back(usedX(fmin(cloc0,cloc2),fmax(cloc0,cloc2)));iu0=used.size()-1; }
259 int iu1 = -1; float w1=0.;
260 for (unsigned int iu=0; iu<used.size(); ++iu) {
261 if (fmin(cloc1,cloc2)==used[iu].locMin && fmax(cloc1,cloc2)==used[iu].locMax) {
262 iu1=iu;
263 if (used[iu].chains.size()>0)
264 for (unsigned int iw=0; iw<used[iu].chains.size(); iw++) w1 += wchain[used[iu].chains[iw].first];
265 //used[iu].add(inext,mch[1]);
266 break;
267 }
268 }
269 if (iu1<0) { used.push_back(usedX(fmin(cloc1,cloc2),fmax(cloc1,cloc2))); iu1=used.size()-1; }
270
271 if ( assoc[cloc2] < 3. && w0 < 2. && w1 < 2.) {
272 accept3[inext] = true;
273 assoc[cloc2] += 1.;
274 used[iu0].add(inext, mchains[inext][1]);
275 used[iu1].add(inext, mchains[inext][1]);
276 if (wchain[inext]==0.5) { // accept the identical chain, too
277 for (unsigned int ic=0; ic< rchains.size(); ++ic) {
278 if (rchains[ic][0]==cloc1 && rchains[ic][1] == cloc0) {
279 accept3[ic]=true;
280 used[iu0].add(ic, mchains[ic][1]);
281 used[iu1].add(ic, mchains[ic][1]);
282 break;
283 }
284 }
285 }
286 }
287 else { // find new association
288 int i1 = rchains[inext][0];
289 int i2 = rchains[inext][1];
290 float mMn = 1000000.;
291 int ipn=-1; int iploc=-1;
292 rchains[inext].pop_back();
293 for (unsigned int i3 = 0; i3 < hs.size(); ++i3) {
294 double q02 = qq(hs[i1],hs[i3],match);
295 if (match > 0.) continue;
296 if (assoc[i3] > 3-wchain[inext]) continue;
297 // check pair association
298 w0=0.; w1=0.;
299 for (unsigned int iu=0; iu<used.size(); ++iu) {
300 if (fmin(cloc0,i3)==used[iu].locMin && fmax(cloc0,i3)==used[iu].locMax ) {
301 if (used[iu].chains.size() > 0)
302 for (unsigned int iw=0; iw<used[iu].chains.size(); ++iw) w0 += wchain[used[iu].chains[iw].first];
303 }
304 if (fmin(cloc1,i3)==used[iu].locMin && fmax(cloc1,i3)==used[iu].locMax ) {
305 if (used[iu].chains.size()>0)
306 for (unsigned int iw=0; iw<used[iu].chains.size(); ++iw) w1 += wchain[used[iu].chains[iw].first];
307 }
308 }
309 if (w0+wchain[inext]>2. || w1+wchain[inext]>2.) continue;
310
311 float q12 = qq(hs[i2],hs[i3],match);
312 float q01 = qq(hs[i1],hs[i2],match);
313 float m = sqrt(9*pim*pim+q02*q02+q01*q01+q12*q12);
314 if (m>0. && m <mMn ) { mMn = m; ipn = i3; iploc = i3; }
315 }
316 if (ipn>=0) {
317 rchains[inext].push_back(ipn); rchains[inext][2]=iploc; mchains[inext][1]=mMn;
318 }
319 else { // chain not recovered
320 wchain[inext]=0.;
321 }
322 }
323 }
324 } // end loop over chains
325 // end 3rd member optimization
326
327 for (unsigned int ip=0; ip < wchain.size(); ++ip) {
328 if (!wchain[ip]) continue;
329 if (rchains[ip].size() < 3) continue;
330 float m3min = mchains[ip][1];
331 if (m3min > 0.59) continue;
332 // dalitz plot
333 std::pair<float,float> dd = dalitz3(hs[rchains[ip][0]], hs[rchains[ip][1]], hs[rchains[ip][2]]);
334 _dalitz->fill(dd.first,dd.second,1.*wchain[ip]);
335 // Delta(Q) spectra
336 float qlmin = mchains[ip][0];
337 float qxmin = qq(hs[rchains[ip][0]], hs[rchains[ip][2]], match);
338 float xlmin = qq(hs[rchains[ip][1]], hs[rchains[ip][2]], match);
339 _Delta3h->fill(qxmin, 0.5*wchain[ip]);
340 _Delta3h->fill(xlmin, 0.5*wchain[ip]);
341 _Delta3h->fill(qlmin, -1.*wchain[ip]);
342 }
343 }
344
345 /// Normalise histograms etc., after the run
346 void finalize() {
347
348
349 // normalize by the number of charged particles
350 // counter automatic division by bin size
351 double norm = 0.01 / (_h_nch->xMean()*_h_nch->numEntries());
352 _dalitz->scaleW(norm);
353 _DeltaQ->scaleW(norm);
354 _Delta3h->scaleW(norm);
355
356 }
357
358 /// @}
359 double qq(const Particle& gp1, const Particle& gp2, bool& match) {
360 match = gp1.charge() * gp2.charge() > 0;
361 FourMomentum p1, p2;
362 p1.setPM(gp1.px(), gp1.py(), gp1.pz(), pim);
363 p2.setPM(gp2.px(), gp2.py(), gp2.pz(), pim);
364 return sqrt(fmax(0., (p1 + p2).mass2() - 4*pim*pim));
365 }
366
367 std::pair<float,float> dalitz3(const Particle& gp1, const Particle& gp2, const Particle& gp3) const {
368
369 float p1= gp1.pt();
370 float p2= gp2.pt();
371 float p3= gp3.pt();
372 float th1 = gp1.theta();
373 float th2 = gp2.theta();
374 float th3 = gp3.theta();
375 float ph1 = gp1.phi();
376 float ph2 = gp2.phi();
377 float ph3 = gp3.phi();
378 float e1 = sqrt(p1*p1+pim*pim);
379 float e2 = sqrt(p2*p2+pim*pim);
380 float e3 = sqrt(p3*p3+pim*pim);
381
382 float p1x = p1*cos(ph1)*sin(th1);
383 float p1y = p1*sin(ph1)*sin(th1);
384 float p1z = p1*cos(th1);
385
386 float p2x = p2*cos(ph2)*sin(th2);
387 float p2y = p2*sin(ph2)*sin(th2);
388 float p2z = p2*cos(th2);
389
390 float p3x = p3*cos(ph3)*sin(th3);
391 float p3y = p3*sin(ph3)*sin(th3);
392 float p3z = p3*cos(th3);
393
394 float px = p1x+p2x+p3x;
395 float py = p1y+p2y+p3y;
396 float pz = p1z+p2z+p3z;
397 float ap = sqrt(px*px+py*py+pz*pz);
398 float e=e1+e2+e3;
399
400 float beta = ap/e;
401 float gamma = 1./sqrt(1-beta*beta);
402
403 float p1l = (p1x*px+p1y*py+p1z*pz)/ap;
404 float p2l = (p2x*px+p2y*py+p2z*pz)/ap;
405 float p3l = (p3x*px+p3y*py+p3z*pz)/ap;
406
407 float e1_boost = gamma*e1-gamma*beta*p1l;
408 float e2_boost = gamma*e2-gamma*beta*p2l;
409 float e3_boost = gamma*e3-gamma*beta*p3l;
410
411 float Q = sqrt(e*e-ap*ap)-3*pim;
412
413 return std::pair<float,float>(sqrt(3.)*(e1_boost-e2_boost)/Q , 3*(e3_boost-pim)/Q-1.);
414 }
415
416 private:
417
418 // Data members like post-cuts event weight counters go here
419 float pim;
420
421 private:
422
423 /// @name Histograms
424 Histo1DPtr _DeltaQ;
425 Histo1DPtr _Delta3h;
426 Histo1DPtr _h_nch;
427 Histo2DPtr _dalitz;
428 /// @}
429 };
430
431 RIVET_DECLARE_PLUGIN(ATLAS_2017_I1624693);
432}
|