Rivet analyses referenceMC_TAUPOLMonte Carlo validation observables for tau polarisationExperiment: () Status: VALIDATED Authors:
Beams: * * Beam energies: ANY Run details:
Observables sensitive to the polarisation of the $\tau$ lepton. Source code: MC_TAUPOL.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/TauFinder.hh"
5
6namespace Rivet {
7
8
9 /// @brief Monte Carlo validation observables for tau polarisation
10 class MC_TAUPOL : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(MC_TAUPOL);
15
16 /// @name Analysis methods
17 //@{
18
19 /// Book histograms and initialise projections before the run
20 void init() {
21
22 _mode = 0; // default is ditau
23 if ( getOption("MODE") == "DITAU" ) { _mode = 0; }
24 if ( getOption("MODE") == "SINGLE" ) { _mode = 1; }
25 _lothresh = getOption("LOMASS", _mode? 75*GeV : 86.*GeV);
26 _hithresh = getOption("HIMASS", _mode? 85*GeV : 96*GeV);
27
28 // Initialise and register projections
29 declare(TauFinder(TauDecay::ANY), "taus");
30
31 // Book histograms
32 // tau invmass
33 book(_h["m_tau"], "h_m_tau", 50,1.,3.);
34 book(_h["m_tautau"], "h_m_tautau", 50,35.,420.);
35 book(_h["m_B"], "h_m_B", 50,35.,420.);
36 book(_h["B_pT"], "h_B_pT", 100,0.,500.);
37
38 // longnitudal spin correlations
39 book(_h["m_tau_vis"], "h_m_tau_vis",50,0.,5.);
40 book(_h["m_tautau_vis"], "h_m_tautau_vis",50,0.,100.);
41 book(_h["m_pipi_vis"], "h_m_pipi_vis",60,0.,100.);
42 book(_h["m_pipi_vis_lm"], "h_m_pipi_vis_lm",60,0.,100.);
43 book(_h["m_pipi_vis_hm"], "h_m_pipi_vis_hm",60,0.,100.);
44
45 // transverse spin correlations
46 book(_h["acoplanarity_plus"], "h_acoplanarity_plus",30,0.,6.2);
47 book(_h["acop_pt40"], "h_acop_pt40",30,0.,6.2);
48 book(_h["acop_pt100"], "h_acop_pt100",30,0.,6.2);
49 book(_h["acop_pt200"], "h_acop_pt200",30,0.,6.2);
50 book(_h["acop_pt1000"], "h_acop_pt1000",30,0.,6.2);
51 book(_h["acom_pt40"], "h_acom_pt40",30,0.,6.2);
52 book(_h["acom_pt100"], "h_acom_pt100",30,0.,6.2);
53 book(_h["acom_pt200"], "h_acom_pt200",30,0.,6.2);
54 book(_h["acom_pt1000"], "h_acom_pt1000",30,0.,6.2);
55 book(_h["acoplanarity_minus"], "h_acoplanarity_minus",30,0.,6.2);
56
57 // leptonic
58
59 book(_h["2B_xlep"], "h_2B_xlep", 50,0.0,1.0);
60 book(_h["2B_xel"], "h_2B_xel", 50,0.0,1.0);
61 book(_h["2B_xmu"], "h_2B_xmu", 50,0.0,1.0);
62 book(_h["xel_B"], "h_xel_B", 50,0.0,1.0);
63 book(_h["xmu_B"], "h_xmu_B", 50,0.0,1.0);
64
65 // hadronic
66
67 book(_h["1B_xpi"], "h_1B_xpi", 30,0.0,1.0);
68 book(_h["1B_xpi_lm"], "h_1B_xpi_lm", 30,0.0,1.0);
69 book(_h["1B_xpi_hm"], "h_1B_xpi_hm", 30,0.0,1.0);
70 book(_h["xpi_B"], "h_xpi_B", 30,0.0,1.0);
71 book(_h["2B_xrho"], "h_2B_xrho", 50,0.0,1.0);
72
73 // angles
74 book(_h["t_pi"], "h_t_pi", 40, -1.0, 1.0);
75 book(_h["t_rho"], "h_t_rho", 40, -1.0, 1.0);
76 book(_h["pipi_theta"], "h_pipi_theta", 40, -1.0, 1.0);
77 book(_ntaus, "h_n_taus", {0, 1, 2, 3, 4});
78
79 }
80
81 void findDecayProducts(const Particle& mother, size_t& nstable,
82 Particles& ep , Particles& em , Particles& nu_e , Particles& nu_ebar,
83 Particles& mup, Particles& mum, Particles& nu_mu, Particles& nu_mubar,
84 Particles& pip, Particles& pim, Particles& pi0,
85 Particles& Kp , Particles& Km , Particles& K0S, Particles& K0L,
86 Particles& eta, Particles& gamma) {
87 for (const Particle& p : mother.children()) {
88 int id = p.pid();
89 if (id == PID::KPLUS) {
90 Kp.push_back(p);
91 ++nstable;
92 }
93 else if (id == PID::KMINUS ) {
94 Km.push_back(p);
95 ++nstable;
96 }
97 else if (id == PID::PIPLUS) {
98 pip.push_back(p);
99 ++nstable;
100 }
101 else if (id == PID::PIMINUS) {
102 pim.push_back(p);
103 ++nstable;
104 }
105 else if (id == PID::EPLUS) {
106 ep.push_back(p);
107 ++nstable;
108 }
109 else if (id == PID::EMINUS) {
110 em.push_back(p);
111 ++nstable;
112 }
113 else if (id == PID::NU_E) {
114 nu_e.push_back(p);
115 ++nstable;
116 }
117 else if (id == PID::NU_EBAR) {
118 nu_ebar.push_back(p);
119 ++nstable;
120 }
121 else if (id == PID::NU_MU) {
122 nu_mu.push_back(p);
123 ++nstable;
124 }
125 else if (id == PID::NU_MUBAR) {
126 nu_mubar.push_back(p);
127 ++nstable;
128 }
129 else if (id == PID::ANTIMUON) {
130 mup.push_back(p);
131 ++nstable;
132 }
133 else if (id == PID::MUON) {
134 mum.push_back(p);
135 ++nstable;
136 }
137 else if (id == PID::PI0) {
138 pi0.push_back(p);
139 ++nstable;
140 }
141 else if (id == PID::K0S) {
142 K0S.push_back(p);
143 ++nstable;
144 }
145 else if (id == PID::K0L) {
146 K0L.push_back(p);
147 ++nstable;
148 }
149 else if (id == PID::ETA) {
150 eta.push_back(p);
151 ++nstable;
152 }
153 else if (id == PID::PHOTON) {
154 gamma.push_back(p);
155 ++nstable;
156 }
157 else if ( !p.children().empty() ) {
158 findDecayProducts(p, nstable, ep, em, nu_e, nu_ebar, mup, mum, nu_mu, nu_mubar,
159 pip, pim, pi0, Kp, Km, K0S, K0L, eta, gamma);
160 }
161 else {
162 ++nstable;
163 }
164 }
165 }
166
167
168 double angleRF(const FourMomentum& p1, const FourMomentum& p2, const FourMomentum& RF) const {
169 LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(RF.betaVec());
170 return cos(boost(p1).angle(boost(p2)));
171 }
172
173 /// Perform the per-event analysis
174 void analyze(const Event& event) {
175
176 Particles taus, t_nu, pions;
177 for (const auto& tau : apply<TauFinder>(event, "taus").particlesByPt()) {
178 //remove any from hadronic decay (only want prompt taus)
179 if (tau.parents()[0].isPrompt()) {
180 taus += tau;
181 _h["m_tau"]->fill(tau.mom().mass()/GeV);
182 }
183 }
184 //veto ditau events without exactly 2 prompt taus
185 if (taus.size() != 2 && _mode == 0) vetoEvent;
186 //veto single-tau events without exactly 1 prompt tau
187 if (taus.size() != 1 && _mode == 1) vetoEvent;
188
189 FourMomentum tautaumom, taunumom, phomom, Bmom;
190 if (taus.size() == 2 && _mode == 0) {
191 tautaumom = taus[0].mom() + taus[1].mom();
192 phomom = sum(scanAncestors(taus,22), Kin::p4, FourMomentum());
193 Bmom = tautaumom + phomom;
194 _h["m_tautau"]->fill(tautaumom.mass());
195 }
196 if (taus.size()==1 && _mode == 1) {
197 t_nu = scanAncestors(taus,16);
198 phomom = sum(scanAncestors(taus,22), Kin::p4, FourMomentum());
199 if(t_nu.size() != 1)vetoEvent;
200 if(taus[0].pid()*t_nu[0].pid() > 0) vetoEvent;
201 taunumom = taus[0].mom() + t_nu[0].mom();
202 Bmom = taunumom + phomom;
203 }
204 _ntaus->fill(int(taus.size()));
205 _h["m_B"]->fill(Bmom.mass());
206 _h["B_pT"]->fill(Bmom.pT());
207 // Check opposite charges
208
209 // Calculate invariant mass of tautau
210 FourMomentum tt_mom_vis;
211 for (const auto& t : taus) {
212 FourMomentum t_mom_vis;
213 for (const auto& c : t.children()) {
214 if (c.isVisible()) {
215 tt_mom_vis += c.mom();
216 t_mom_vis += c.mom();
217 }
218 }
219 _h["m_tau_vis"]->fill(t_mom_vis.mass());
220 }
221
222 _h["m_tautau_vis"]->fill(tt_mom_vis.mass());
223 size_t OS = 2;
224 if (Bmom.mass() < _lothresh) OS = 0;
225 else if (Bmom.mass() < _hithresh) OS = 1;
226
227 for (const Particle& tau :taus) {
228 size_t nstable(0);
229 Particles ep,em,nu_e,nu_ebar,mup,mum,nu_mu,nu_mubar;
230 Particles pip, pim, pi0, Kp , Km, K0S, K0L, eta,gamma;
231 findDecayProducts(tau, nstable,ep,em,nu_e,nu_ebar,mup,mum,nu_mu,nu_mubar,
232 pip, pim, pi0, Kp, Km, K0S, K0L,eta,gamma);
233 if (tau.pid() < 0) {
234 swap(pim,pip);
235 swap(Kp,Km);
236 swap(em,ep);
237 swap(mum,mup);
238 swap(nu_e ,nu_ebar );
239 swap(nu_mu,nu_mubar);
240 }
241 // 2 hadrons
242 FourMomentum tmom = tau.mom();
243 const LorentzTransform tboost = LorentzTransform::mkFrameTransformFromBeta(tmom.betaVec());
244 if (nstable==2) {
245 if (pim.size()==1) {
246 const double xpi = pim[0].mom().p()/tau.mom().p();
247 const FourMomentum pimom = pim[0].mom();
248 const double angle = cos(tboost(pimom).angle(tmom));
249 _h["xpi_B"]->fill(pim[0].mom().p()/Bmom.p());
250 _h["t_pi"]->fill(angle);
251 if (OS==0) _h["1B_xpi_lm"]->fill(xpi);
252 else if (OS==1) _h["1B_xpi"]->fill(xpi);
253 else _h["1B_xpi_hm"]->fill(xpi);
254 pions.push_back(pim[0]);
255 }
256 }
257 else if (nstable==3 ) {
258 if (em.size()==1 && nu_ebar.size()==1) {
259 const double xlep = em[0].mom().p()/tau.mom().p();
260 _h["2B_xlep"]->fill(xlep);
261 _h["2B_xel"]->fill(xlep);
262 _h["xel_B"]->fill(em[0].mom().p()/Bmom.p());
263 }
264 else if (mum.size()==1 && nu_mubar.size()==1) {
265 const double xlep = mum[0].mom().p()/tau.mom().p();
266 _h["2B_xlep"]->fill(xlep);
267 _h["2B_xmu"]->fill(xlep);
268 _h["xmu_B"]->fill(mum[0].mom().p()/Bmom.p());
269 }
270 else if (pim.size()==1 && pi0.size()==1) {
271 const FourMomentum ptot = pim[0].mom()+pi0[0].mom();
272 const double xrho = ptot.p()/tau.mom().p();
273 _h["2B_xrho"]->fill(xrho);
274 const double angle = cos(tboost(ptot).angle(tmom));
275 _h["t_rho"]->fill(angle);
276 }
277 } //end nstable = 3
278
279 } //end tau loop
280
281 // TAU LONGNITUDAL SPIN EFFECTS PART
282 //fill tautau->pipi inv mass
283 if (pions.size() == 2){
284 FourMomentum pi1mom = pions[0].mom();
285 FourMomentum pi2mom = pions[1].mom();
286 FourMomentum pipi = pions[0].mom() + pions[1].mom();
287 if (OS==0) _h["m_pipi_vis_lm"]->fill(pipi.mass()/GeV);
288 else if (OS==1) _h["m_pipi_vis"]->fill(pipi.mass()/GeV);
289 else _h["m_pipi_vis_hm"]->fill(pipi.mass()/GeV);
290 _h["pipi_theta"]->fill(angleRF(pi1mom,pi2mom,Bmom));
291 }
292
293 // TAU TRANSVERSE SPIN EFFECTS PART
294 if (pions.size() == 2 && _mode != 1 && OS == 1){
295
296 FourMomentum tautau = taus[0].mom() + taus[1].mom();
297
298 // calculate angles between (tau, tau_pi) and (tau2, tau2_pi) planes
299 const LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(tautau.betaVec());
300 taus[0] = Particle(PID::ANY, boost.transform(taus[0].mom()));
301 pions[0] = Particle(PID::ANY, boost.transform(pions[0].mom()));
302 taus[1] = Particle(PID::ANY, boost.transform(taus[1].mom()));
303 pions[1] = Particle(PID::ANY, boost.transform(pions[1].mom()));
304
305 // calculate the angle
306 Vector3 t0xp0 = taus[0].p3().cross(pions[0].p3()).unit();
307 Vector3 t1xp1 = taus[1].p3().cross(pions[1].p3()).unit();
308 double phi_star_sign = acos(t0xp0.dot(t1xp1));
309
310 // extend acoplanarity angle to range (0, 2pi)
311 Vector3 p0xp1 = pions[0].p3().cross(pions[1].p3()).unit();
312 double sign = p0xp1.dot(taus[0].p3().unit());
313 if (sign > 0) phi_star_sign = 2* M_PI - phi_star_sign;
314
315
316 // alpha angles in the lab frame
317 Vector3 tlv3_Ez(0.0, 0.0, 1.0 );
318 Vector3 cross1_pi = tlv3_Ez.cross(taus[0].p3()).unit();
319 Vector3 cross2_pi = pions[0].p3().cross(taus[0].p3()).unit();
320 const double alpha_pi = acos(abs(cross1_pi.dot(cross2_pi)));
321
322 if (alpha_pi > M_PI/4.) {
323 _h["acoplanarity_plus"]->fill(phi_star_sign);
324 if (Bmom.pT() < 40*GeV) _h["acop_pt40"]->fill(phi_star_sign);
325 else if (Bmom.pT() < 100*GeV) _h["acop_pt100"]->fill(phi_star_sign);
326 else if (Bmom.pT() < 200*GeV) _h["acop_pt200"]->fill(phi_star_sign);
327 else _h["acop_pt1000"]->fill(phi_star_sign);
328 }
329
330 if (alpha_pi < M_PI/4.) {
331 _h["acoplanarity_minus"]->fill(phi_star_sign);
332 if (Bmom.pT() < 40*GeV) _h["acom_pt40"]->fill(phi_star_sign);
333 else if (Bmom.pT() < 100*GeV) _h["acom_pt100"]->fill(phi_star_sign);
334 else if (Bmom.pT() < 200*GeV) _h["acom_pt200"]->fill(phi_star_sign);
335 else _h["acom_pt1000"]->fill(phi_star_sign);
336 }
337
338 } // end of tau spin corr
339
340 } //end void function
341
342
343 /// Normalise histograms etc., after the run
344 void finalize() {
345 normalize(_h);
346 normalize(_ntaus);
347 }
348
349 //@}
350
351 Particles scanAncestors(const Particles& taus, int id){
352 Particles rtn;
353 for (const Particle& tau : taus) {
354 Particle thistau = tau;
355 while (true) {
356 // check the tau has parents first else end loop
357 if (thistau.parents().empty()) break;
358 //loop over taus siblings
359 for (const Particle& sibling : thistau.parents()[0].children()) {
360 // find ID
361 if (sibling.abspid() == id) {
362 // add particles to rtn vector
363 bool duplicate = true;
364 for (const auto& p : rtn) {
365 if (sibling.isSame(p)) duplicate = false;
366 }
367 if (duplicate) rtn += sibling;
368 }
369 }
370 // want to loop over taus parent's siblings next
371 thistau = thistau.parents()[0];
372 // if taus parent is not a tau end the loop
373 if (thistau.abspid() != PID::TAU) break;
374 }
375 }
376 return rtn;
377 }
378
379 /// @name Histograms
380 //@{
381
382 // histograms for leptonic decay
383 map<string,Histo1DPtr> _h;
384
385 BinnedHistoPtr<int> _ntaus;
386
387 size_t _mode;
388
389 double _lothresh, _hithresh;
390
391 //@}
392
393 };
394
395
396 RIVET_DECLARE_PLUGIN(MC_TAUPOL);
397}
|