Rivet analyses referenceDELPHI_2000_I522656Events shapes at MZ as function of thrust directionExperiment: DELPHI (LEP) Inspire ID: 522656 Status: VALIDATED Authors:
Beam energies: (45.6, 45.6) GeV Run details:
Measurement of event shapes for different ranges of the angle between the thrust axis and the beam. Jet rate using a mass measure or the Geneva algorithm are not implemented. Source code: DELPHI_2000_I522656.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/FinalState.hh"
5#include "Rivet/Projections/FastJets.hh"
6#include "Rivet/Projections/Thrust.hh"
7#include "Rivet/Projections/Sphericity.hh"
8#include "Rivet/Projections/Hemispheres.hh"
9#include "Rivet/Projections/ParisiTensor.hh"
10#include "Rivet/Tools/BinnedHistogram.hh"
11#include "fastjet/EECambridgePlugin.hh"
12
13namespace Rivet {
14
15
16 /// @brief event shapes vs thrust direction
17 class DELPHI_2000_I522656 : public Analysis {
18 public:
19
20 /// Constructor
21 RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_2000_I522656);
22
23
24 /// @name Analysis methods
25 ///@{
26
27 /// Book histograms and initialise projections before the run
28 void init() {
29 // Initialise and register projections.
30 declare(Beam(), "Beams");
31 const FinalState fs;
32 declare(fs, "FS");
33 const Thrust thrust(fs);
34 declare(thrust, "Thrust");
35 declare(Sphericity(fs), "Sphericity");
36 declare(ParisiTensor(fs), "Parisi");
37 declare(Hemispheres(thrust), "Hemispheres");
38 declare(FastJets(fs, FastJets::DURHAM, 0.7), "DurhamJets");
39 declare(FastJets(fs, FastJets::JADE , 0.7), "JadeJets" );
40
41 // book histograms
42 vector<double> bins={0.00,0.12,0.24,0.36,0.48,0.60,0.72,0.84,0.96};
43 // thrust angle binned
44 unsigned int iy=1,ioff=0;
45 for(unsigned int ix=0;ix<8;++ix) {
46 {Histo1DPtr tmp; _h_EEC.add(bins[ix],bins[ix+1],book(tmp,21+ioff,1,iy));}
47 {Histo1DPtr tmp;_h_AEEC.add(bins[ix],bins[ix+1],book(tmp,25+ioff,1,iy));}
48 {Histo1DPtr tmp;_h_cone.add(bins[ix],bins[ix+1],book(tmp,29+ioff,1,iy));}
49 if(ioff==0) {
50 Histo1DPtr tmp;_h_thrust.add(bins[iy],bins[iy+1],book(tmp,33,1,iy));
51 }
52 ++iy;
53 if(iy==3) {
54 ++ioff;
55 iy=1;
56 }
57 }
58 // total values
59 book(_h_EEC_all , 3,1,1);
60 book(_h_AEEC_all , 4,1,1);
61 book(_h_cone_all , 5,1,1);
62 book(_h_thrust_all, 6,1,1);
63 book(_h_Oblateness, 7,1,1);
64 book(_h_C , 8,1,1);
65 book(_h_heavy , 9,1,1);
66 book(_h_sum ,10,1,1);
67 book(_h_diff ,11,1,1);
68 book(_h_wide ,12,1,1);
69 book(_h_total ,13,1,1);
70 book(_h_jade ,17,1,1);
71 book(_h_dur ,18,1,1);
72 book(_h_cam ,20,1,1);
73 book(_h_bin,"/TMP/hbin",bins);
74 }
75
76
77 /// Perform the per-event analysis
78 void analyze(const Event& event) {
79 const FinalState& fs = apply<FinalState>(event, "FS");
80 if ( fs.particles().size() < 2) vetoEvent;
81 // Get beams and average beam momentum
82 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
83 // thrust
84 const Thrust& thrust = apply<Thrust>(event, "Thrust");
85 // angle bettwen thrust and beam
86 double cosThrust = abs(beams.first.p3().unit().dot(thrust.thrustAxis()));
87 _h_bin->fill(cosThrust);
88 // thrust and related
89 _h_thrust_all->fill( 1.-thrust.thrust());
90 _h_thrust .fill(cosThrust, 1.-thrust.thrust());
91 _h_Oblateness->fill(thrust.oblateness() );
92
93 // visible energy and make pseudojets
94 double Evis = 0.0;
95 PseudoJets pjs;
96 for (const Particle& p : fs.particles()) {
97 Evis += p.E();
98 fastjet::PseudoJet pj = p;
99 pjs.push_back(pj);
100 }
101 double Evis2 = sqr(Evis);
102 // (A)EEC
103 // Need iterators since second loop starts at current outer loop iterator, i.e. no "foreach" here!
104 for (Particles::const_iterator p_i = fs.particles().begin(); p_i != fs.particles().end(); ++p_i) {
105 for (Particles::const_iterator p_j = p_i; p_j != fs.particles().end(); ++p_j) {
106 if(p_i == p_j) continue;
107 const Vector3 mom3_i = p_i->momentum().p3();
108 const Vector3 mom3_j = p_j->momentum().p3();
109 const double energy_i = p_i->momentum().E();
110 const double energy_j = p_j->momentum().E();
111 const double thetaij = 180.*mom3_i.unit().angle(mom3_j.unit())/M_PI;
112 double eec = (energy_i*energy_j) / Evis2;
113 eec *= 2.;
114 _h_EEC_all->fill( thetaij, eec);
115 _h_EEC .fill(cosThrust, thetaij, eec);
116 if (thetaij <90.){
117 _h_AEEC_all->fill( thetaij, -eec);
118 _h_AEEC .fill(cosThrust, thetaij, -eec);
119 }
120 else {
121 _h_AEEC_all->fill( 180.-thetaij, eec);
122 _h_AEEC .fill(cosThrust,180.-thetaij, eec);
123 }
124 }
125 }
126 // hemisphere related
127 const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
128 _h_heavy->fill(hemi.scaledM2high());
129 _h_diff ->fill(hemi.scaledM2diff());
130 _h_sum ->fill(hemi.scaledM2low()+hemi.scaledM2high());
131 _h_wide ->fill(hemi.Bmax() );
132 _h_total->fill(hemi.Bsum() );
133 // C-parameter
134 const ParisiTensor& parisi = apply<ParisiTensor>(event, "Parisi");
135 _h_C->fill(parisi.C());
136 // jets
137 const FastJets& durjet = apply<FastJets>(event, "DurhamJets");
138 const FastJets& jadejet = apply<FastJets>(event, "JadeJets");
139 if (durjet .clusterSeq()) _h_dur ->fill( durjet.clusterSeq()->exclusive_ymerge_max(2));
140 if (jadejet.clusterSeq()) _h_jade->fill(jadejet.clusterSeq()->exclusive_ymerge_max(2));
141 // Cambridge is more complicated, inclusive defn
142 for (size_t i = 0; i < _h_cam->numBins(); ++i) {
143 double ycut = _h_cam->bin(i).xMax();
144 // double width = _h_y_2_Cambridge->bin(i).width();
145 fastjet::EECambridgePlugin plugin(ycut);
146 fastjet::JetDefinition jdef(&plugin);
147 fastjet::ClusterSequence cseq(pjs, jdef);
148 unsigned int njet = cseq.inclusive_jets().size();
149 if(njet==2) {
150 _h_cam->fill( _h_cam->bin(i).xMid());
151 break;
152 }
153 }
154 // jet cone
155 Vector3 jetAxis=thrust.thrustAxis();
156 if(hemi.highMassDirection()) jetAxis *=-1.;
157 for(const Particle & p : fs.particles()) {
158 const double thetaij = 180.*jetAxis.angle(p.p3().unit())/M_PI;
159 double jcef = p.E()/ Evis;
160 _h_cone_all->fill( thetaij,jcef);
161 _h_cone .fill(cosThrust,thetaij,jcef);
162 }
163 }
164
165
166 /// Normalise histograms etc., after the run
167 void finalize() {
168 for(unsigned int ix=0;ix<8;++ix) {
169 if(ix<2) scale(_h_thrust.histos()[ix],1./_h_bin->bins()[ix].area());
170 scale(_h_EEC.histos()[ix],180./M_PI/_h_bin->bins()[ix].area());
171 scale(_h_AEEC.histos()[ix],180./M_PI/_h_bin->bins()[ix].area());
172 scale(_h_cone.histos()[ix],180./M_PI/_h_bin->bins()[ix].area());
173 }
174
175
176 // _h_thrust.scale(1./sumOfWeights(),this);
177 // _h_EEC.scale(180./M_PI/sumOfWeights(),this);
178 // _h_AEEC.scale(180./M_PI/sumOfWeights(),this);
179 // _h_cone.scale(180./M_PI/sumOfWeights(),this);
180
181
182
183
184 scale(_h_thrust_all, 1./sumOfWeights());
185 scale(_h_EEC_all, 180./M_PI/sumOfWeights());
186 scale(_h_AEEC_all, 180./M_PI/sumOfWeights());
187 scale(_h_cone_all, 180./M_PI/sumOfWeights());
188 scale(_h_Oblateness, 1./sumOfWeights());
189 scale(_h_C , 1./sumOfWeights());
190 scale(_h_heavy , 1./sumOfWeights());
191 scale(_h_sum , 1./sumOfWeights());
192 scale(_h_diff , 1./sumOfWeights());
193 scale(_h_wide , 1./sumOfWeights());
194 scale(_h_total , 1./sumOfWeights());
195 scale(_h_dur , 1./sumOfWeights());
196 scale(_h_jade , 1./sumOfWeights());
197 scale(_h_cam , 1./sumOfWeights());
198 }
199
200 ///@}
201
202
203 /// @name Histograms
204 ///@{
205 Histo1DPtr _h_thrust_all,_h_EEC_all,_h_AEEC_all,_h_cone_all;
206 Histo1DPtr _h_Oblateness,_h_C,_h_heavy,_h_sum,_h_diff,_h_wide,_h_total;
207 Histo1DPtr _h_jade,_h_dur,_h_cam;
208 BinnedHistogram _h_thrust,_h_EEC,_h_AEEC,_h_cone;
209 Histo1DPtr _h_bin;
210 ///@}
211
212
213 };
214
215
216 RIVET_DECLARE_PLUGIN(DELPHI_2000_I522656);
217
218}
|