Rivet analyses referenceATLAS_2011_I914491Long-lived heavy charged particle searchExperiment: ATLAS (LHC) Inspire ID: 914491 Status: UNVALIDATED Authors:
Beam energies: (3500.0, 3500.0) GeV Run details:
ATLAS search for long-lived heavy charged particles for four different mass cuts. Currently only the slepton search is implemented. Source code: ATLAS_2011_I914491.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Math/Constants.hh"
4#include "Rivet/Projections/ChargedFinalState.hh"
5#include "Rivet/Projections/NonHadronicFinalState.hh"
6#include "Rivet/Projections/VetoedFinalState.hh"
7
8namespace Rivet {
9
10
11 class ATLAS_2011_I914491 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_I914491);
16
17
18 /// @name Analysis methods
19 /// @{
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23
24 // get the non-hadronic final-state particles
25 double etaMax = 2.5;
26 const NonHadronicFinalState nhfs((Cuts::etaIn(-etaMax,etaMax) && Cuts::pT >= 13.*GeV));
27 declare(nhfs,"NHFS");
28 // select the charged ones
29 const ChargedFinalState cfs(nhfs);
30 declare(cfs,"CFS");
31 // and then veto electrons, and taus to be safe
32 VetoedFinalState vfs(cfs);
33 vfs.addVetoPairId(PID::ELECTRON);
34
35 declare(vfs,"VFS");
36
37 /// Book histograms
38 book(_count_trigger ,"count_trigger" , 1, 0., 1.);
39 book(_count_event ,"count_selection", 1, 0., 1.);
40 book(_count_quality ,"count_quality" , 1, 0., 1.);
41 book(_count_beta ,"count_beta" , 1, 0., 1.);
42 book(_count_90 ,"count_90" , 1, 0., 1.);
43 book(_count_110 ,"count_110", 1, 0., 1.);
44 book(_count_120 ,"count_120", 1, 0., 1.);
45 book(_count_130 ,"count_130", 1, 0., 1.);
46
47 book(_hist_beta ,"beta",1000, 0., 2.);
48 book(_hist_time ,"time",1000, -50, 50.);
49 book(_hist_mass ,"mass", 60, 5., 305.);
50 }
51
52
53 double rndGauss(double sigma, double mean) {
54 double r = sqrt(-2.0*log(rand()/static_cast<double>(RAND_MAX)));
55 double phi = rand()/static_cast<double>(RAND_MAX)*2.0*pi;
56 return mean+sigma*r*sin(phi);
57 }
58
59 /// Perform the per-event analysis
60 void analyze(const Event& event) {
61 // smearing parameters
62 // time measurement (in ns)
63// const double tsmear=5.*0.7;
64 const double tsmear=0.7;
65 // sagita error
66 const double csag =1.1E-4;
67 // multiple scattering
68 const double cms =2.0E-2;
69 // muon chamber radius (in metres)
70 const double radius = 10.e3;
71 // convert to ns
72 const double tr = radius/c_light;
73 // get the charged final-state particles
74 Particles charged =
75 apply<VetoedFinalState>(event,"VFS").particles();
76 // need at least two candidates
77 if(charged.size()<2) vetoEvent;
78 // number passing trigger
79 _count_trigger->fill(0.5);
80 // Z mass veto
81 for ( const Particle & mu1 : charged ) {
82 for ( const Particle & mu2 : charged ) {
83 double mass = (mu1.momentum()+mu2.momentum()).mass();
84 double diff = abs(mass-91.18);
85 if(diff<10.) vetoEvent;
86 }
87 }
88 // number passing first event selection
89 _count_event->fill(0.5);
90 // now find the candidates
91 // loop over the particles and find muons and heavy charged particles
92 map<double,Particle> muonCandidates;
93 for (const Particle& mu : charged) {
94 // calculate the smeared momentum
95 double pT = mu.pT2();
96 double pmag = sqrt(pT+sqr(mu.pz()));
97 double deltap = sqrt( sqr(csag*sqr(pmag)) + sqr(cms*mu.E()/GeV));
98 double psmear = rndGauss(deltap,pmag);
99 // keep particles with pT>40
100 if(psmear/pmag*mu.perp()<40.*GeV||
101 psmear/pmag*mu.perp()>1000.*GeV) continue;
102 muonCandidates.insert(make_pair(psmear,mu));
103 }
104 // require two candidates
105 if(muonCandidates.size()<2) vetoEvent;
106 // number passing "quality" cut
107 _count_quality->fill(0.5);
108 // now do the time of flight
109 bool filled = false;
110 for(map<double,Particle>::const_iterator it=muonCandidates.begin();
111 it!=muonCandidates.end();++it) {
112 // true magnitude and pT of momentum
113 double pT = it->second.pT2();
114 double pmag = sqrt(pT+sqr(it->second.pz()));
115 pT = sqrt(pT);
116 // true time difference in ns
117 double deltaT =tr *(it->second.E()-pmag)/pT;
118 // smear it
119 deltaT = rndGauss(tsmear,deltaT);
120 // beta
121 double beta = 1./(1.+deltaT/tr*pT/pmag);
122 _hist_beta->fill(beta);
123 _hist_time->fill(deltaT);
124 // beta cut
125 if(beta<0.95) continue;
126 // mass
127 double mass = 2.*pT*it->first*deltaT/tr*(1.+0.5*deltaT/tr*pT/pmag);
128 if(mass<0.) continue;
129 mass = sqrt(mass);
130 filled = true;
131 _hist_mass->fill(mass);
132 if(mass>90. ) {
133 _count_90 ->fill(0.5);
134 if(mass>110.) {
135 _count_110->fill(0.5);
136 if(mass>120.) {
137 _count_120->fill(0.5);
138 if(mass>130.) {
139 _count_130->fill(0.5);
140 }
141 }
142 }
143 }
144 }
145 if(!filled) vetoEvent;
146 // number passing beta cut
147 _count_beta->fill(0.5);
148 }
149
150 /// @}
151
152 void finalize() {
153 double fact = crossSection()/picobarn/sumOfWeights()*37;
154 MSG_WARNING("testing " << crossSection()/picobarn << " " << sumOfWeights() << " " << fact);
155 scale(_hist_beta,fact);
156 scale(_hist_time,fact);
157 scale(_hist_mass,fact);
158 scale( _count_trigger, fact);
159 scale( _count_event, fact);
160 scale( _count_quality, fact);
161 scale( _count_beta, fact);
162 scale( _count_90, fact);
163 scale( _count_110, fact);
164 scale( _count_120, fact);
165 scale( _count_130, fact);
166 }
167
168 private:
169
170 /// @name Histograms
171 /// @{
172 Histo1DPtr _hist_beta;
173 Histo1DPtr _hist_time;
174 Histo1DPtr _hist_mass;
175 Histo1DPtr _count_trigger;
176 Histo1DPtr _count_event;
177 Histo1DPtr _count_quality;
178 Histo1DPtr _count_beta;
179 Histo1DPtr _count_90;
180 Histo1DPtr _count_110;
181 Histo1DPtr _count_120;
182 Histo1DPtr _count_130;
183 /// @}
184
185 };
186
187
188 RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_I914491, ATLAS_2011_S9108483);
189
190}
|