1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/IdentifiedFinalState.hh"
  5#include "Rivet/Projections/PromptFinalState.hh"
  6#include "Rivet/Projections/LeptonFinder.hh"
  7#include "Rivet/Projections/PartonicTops.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// CMS 8 TeV dilepton channel ttbar spin correlations and polarisation analysis
 13  class CMS_2016_I1413748 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2016_I1413748);
 18
 19
 20    /// Book histograms and initialise projections
 21    void init() {
 22
 23      // Complete final state
 24      FinalState fs;
 25
 26      // Projection for dressed electrons and muons
 27      IdentifiedFinalState photons(fs);
 28      photons.acceptIdPair(PID::PHOTON);
 29
 30      IdentifiedFinalState el_id(fs);
 31      el_id.acceptIdPair(PID::ELECTRON);
 32      PromptFinalState electrons(el_id);
 33      declare(electrons, "Electrons");
 34      LeptonFinder dressed_electrons(electrons, photons, 0.1);
 35      declare(dressed_electrons, "DressedElectrons");
 36
 37      IdentifiedFinalState mu_id(fs);
 38      mu_id.acceptIdPair(PID::MUON);
 39      PromptFinalState muons(mu_id);
 40      declare(muons, "Muons");
 41      LeptonFinder dressed_muons(muons, photons, 0.1);
 42      declare(dressed_muons, "DressedMuons");
 43
 44      // Parton-level top quarks
 45      declare(PartonicTops(TopDecay::E_MU, PromptEMuFromTau::NO), "LeptonicPartonTops");
 46
 47
 48      // Booking of histograms
 49
 50      // This histogram is independent of the parton-level information, and is an addition to the original analysis.
 51      // It is compared to the same data as the parton-level delta_phi histogram d02-x01-y01.
 52      book(_h_dphidressedleptons, "d00-x01-y01", _bins_dphi);
 53
 54      // The remaining histos use parton-level information
 55      book(_h_dphi, "d02-x01-y01", _bins_dphi);
 56      book(_h_cos_opening_angle, "d05-x01-y01", _bins_cos_opening_angle);
 57      book(_h_c1c2, "d08-x01-y01", _bins_c1c2);
 58      book(_h_lep_costheta, "d11-x01-y01", _bins_lep_costheta);
 59      book(_h_lep_costheta_CPV, "d14-x01-y01", _bins_lep_costheta_CPV);
 60
 61      // 2D histos
 62      book(_h_dphi_var[0], "d20-x01-y01", _bins_dphi, _bins_tt_mass);
 63      book(_h_cos_opening_angle_var[0], "d26-x01-y01", _bins_cos_opening_angle, _bins_tt_mass);
 64      book(_h_c1c2_var[0], "d32-x01-y01", _bins_c1c2, _bins_tt_mass);
 65      book(_h_lep_costheta_var[0], "d38-x01-y01", _bins_lep_costheta, _bins_tt_mass);
 66      book(_h_lep_costheta_CPV_var[0], "d44-x01-y01", _bins_lep_costheta_CPV, _bins_tt_mass);
 67
 68      book(_h_dphi_var[1], "d50-x01-y01", _bins_dphi, _bins_tt_pT);
 69      book(_h_cos_opening_angle_var[1], "d56-x01-y01", _bins_cos_opening_angle, _bins_tt_pT);
 70      book(_h_c1c2_var[1], "d62-x01-y01", _bins_c1c2, _bins_tt_pT);
 71      book(_h_lep_costheta_var[1], "d68-x01-y01", _bins_lep_costheta, _bins_tt_pT);
 72      book(_h_lep_costheta_CPV_var[1], "d74-x01-y01", _bins_lep_costheta_CPV, _bins_tt_pT);
 73
 74      book(_h_dphi_var[2], "d80-x01-y01", _bins_dphi, _bins_tt_absrapidity);
 75      book(_h_cos_opening_angle_var[2], "d86-x01-y01", _bins_cos_opening_angle, _bins_tt_absrapidity);
 76      book(_h_c1c2_var[2], "d92-x01-y01", _bins_c1c2, _bins_tt_absrapidity);
 77      book(_h_lep_costheta_var[2], "d98-x01-y01", _bins_lep_costheta, _bins_tt_absrapidity);
 78      book(_h_lep_costheta_CPV_var[2], "d104-x01-y01", _bins_lep_costheta_CPV, _bins_tt_absrapidity);
 79
 80      // Profile histos for asymmetries
 81      book(_h_dphi_profile[0], "d17-x01-y01", _bins_tt_mass);
 82      book(_h_cos_opening_angle_profile[0], "d23-x01-y01", _bins_tt_mass);
 83      book(_h_c1c2_profile[0], "d29-x01-y01", _bins_tt_mass);
 84      book(_h_lep_costheta_profile[0], "d35-x01-y01", _bins_tt_mass);
 85      book(_h_lep_costheta_CPV_profile[0], "d41-x01-y01", _bins_tt_mass);
 86
 87      book(_h_dphi_profile[1], "d47-x01-y01", _bins_tt_pT);
 88      book(_h_cos_opening_angle_profile[1], "d53-x01-y01", _bins_tt_pT);
 89      book(_h_c1c2_profile[1], "d59-x01-y01", _bins_tt_pT);
 90      book(_h_lep_costheta_profile[1], "d65-x01-y01", _bins_tt_pT);
 91      book(_h_lep_costheta_CPV_profile[1], "d71-x01-y01", _bins_tt_pT);
 92
 93      book(_h_dphi_profile[2], "d77-x01-y01", _bins_tt_absrapidity);
 94      book(_h_cos_opening_angle_profile[2], "d83-x01-y01", _bins_tt_absrapidity);
 95      book(_h_c1c2_profile[2], "d89-x01-y01", _bins_tt_absrapidity);
 96      book(_h_lep_costheta_profile[2], "d95-x01-y01", _bins_tt_absrapidity);
 97      book(_h_lep_costheta_CPV_profile[2], "d101-x01-y01", _bins_tt_absrapidity);
 98
 99    }
100
101
102    /// Perform the per-event analysis
103    void analyze(const Event& event) {
104
105
106      // Use particle-level leptons for the first histogram
107      const LeptonFinder& dressed_electrons = apply<LeptonFinder>(event, "DressedElectrons");
108      const LeptonFinder& dressed_muons = apply<LeptonFinder>(event, "DressedMuons");
109
110      const DressedLeptons dressedels = dressed_electrons.dressedLeptons();
111      const DressedLeptons dressedmus = dressed_muons.dressedLeptons();
112
113      const size_t ndressedel = dressedels.size();
114      const size_t ndressedmu = dressedmus.size();
115
116      // For the particle-level histogram, require exactly one electron and exactly one muon, to select
117      // the ttbar->emu channel. Note this means ttbar->emu events with additional PromptFinalState
118      // dilepton pairs from the shower are vetoed - for PYTHIA8, this affects ~0.5% of events, so the
119      // effect is well below the level of sensitivity of the measured distribution.
120      if ( ndressedel == 1 && ndressedmu == 1 ) {
121
122        const int electrontouse = 0, muontouse = 0;
123
124        // Opposite-charge leptons only
125        if ( sameSign(dressedels[electrontouse],dressedmus[muontouse]) ) {
126          MSG_INFO("Error, e and mu have same charge, skipping event");
127        }
128        else {
129          //Get the four-momenta of the positively- and negatively-charged leptons
130          FourMomentum lepPlus = dressedels[electrontouse].charge() > 0 ? dressedels[electrontouse] : dressedmus[muontouse];
131          FourMomentum lepMinus = dressedels[electrontouse].charge() > 0 ? dressedmus[muontouse] : dressedels[electrontouse];
132
133          // Now calculate the variable
134          double dphi_temp = deltaPhi(lepPlus,lepMinus);
135
136          fillWithUFOF( _h_dphidressedleptons, dphi_temp );
137        }
138
139      }
140
141
142      // The remaining variables use parton-level information.
143
144      // Get the leptonically decaying tops
145      const Particles& leptonicpartontops = apply<ParticleFinder>(event, "LeptonicPartonTops").particlesByPt();
146      Particles chargedleptons;
147      unsigned int ntrueleptonictops = 0;
148      bool oppositesign = false;
149
150      if ( leptonicpartontops.size() == 2 ) {
151        for (size_t k = 0; k < leptonicpartontops.size(); ++k) {
152
153          // Get the lepton
154          const Particle lepTop = leptonicpartontops[k];
155          const auto isPromptChargedLepton = [](const Particle& p){return (isChargedLepton(p) && isPrompt(p, false, false));};
156          Particles lepton_candidates = lepTop.allDescendants(firstParticleWith(isPromptChargedLepton), false);
157          if ( lepton_candidates.size() < 1 ) MSG_WARNING("error, TopDecay::E_MU top quark had no daughter lepton candidate, skipping event.");
158
159          // In some cases there is no lepton from the W decay but only leptons from the decay of a radiated gamma.
160          // These hadronic PartonicTops are currently being mistakenly selected by TopDecay::E_MU (as of April 2017), and need to be rejected.
161          // TopDecay::E_MU is being fixed in Rivet, and when it is the veto below should do nothing.
162          /// @todo Should no longer be necessary -- remove
163          bool istrueleptonictop = false;
164          for (size_t i = 0; i < lepton_candidates.size(); ++i) {
165            const Particle& lepton_candidate = lepton_candidates[i];
166            if ( lepton_candidate.hasParentWith(Cuts::pid == PID::PHOTON) ) {
167              MSG_DEBUG("Found gamma parent, top: " << k+1 << " of " << leptonicpartontops.size() << " , lepton: " << i+1 << " of " << lepton_candidates.size());
168              continue;
169            }
170            if ( !istrueleptonictop && sameSign(lepTop,lepton_candidate) ) {
171              chargedleptons.push_back(lepton_candidate);
172              istrueleptonictop = true;
173            }
174            else MSG_WARNING("Found extra prompt charged lepton from top decay (and without gamma parent), ignoring it.");
175          }
176          if ( istrueleptonictop ) ++ntrueleptonictops;
177        }
178      }
179
180      if ( ntrueleptonictops == 2 ) {
181        oppositesign = !( sameSign(chargedleptons[0],chargedleptons[1]) );
182        if ( !oppositesign ) MSG_WARNING("error, same charge tops, skipping event.");
183      }
184
185      if ( ntrueleptonictops == 2 && oppositesign ) {
186
187        // Get the four-momenta of the positively- and negatively-charged leptons
188        FourMomentum lepPlus = chargedleptons[0].charge() > 0 ? chargedleptons[0] : chargedleptons[1];
189        FourMomentum lepMinus = chargedleptons[0].charge() > 0 ? chargedleptons[1] : chargedleptons[0];
190
191        const double dphi_temp = deltaPhi(lepPlus,lepMinus);
192
193        // Get the four-momenta of the positively- and negatively-charged tops
194        FourMomentum topPlus_p4 = leptonicpartontops[0].pid() > 0 ? leptonicpartontops[0] : leptonicpartontops[1];
195        FourMomentum topMinus_p4 = leptonicpartontops[0].pid() > 0 ? leptonicpartontops[1] : leptonicpartontops[0];
196        const FourMomentum ttbar_p4 = topPlus_p4 + topMinus_p4;
197
198        const double tt_mass_temp = ttbar_p4.mass();
199        const double tt_absrapidity_temp = ttbar_p4.absrapidity();
200        const double tt_pT_temp = ttbar_p4.pT();
201
202        // Lorentz transformations to calculate the spin observables in the helicity basis
203
204        // Transform everything to the ttbar CM frame
205        LorentzTransform ttCM;
206        ttCM.setBetaVec(-ttbar_p4.betaVec());
207
208        topPlus_p4 = ttCM.transform(topPlus_p4);
209        topMinus_p4 = ttCM.transform(topMinus_p4);
210
211        lepPlus = ttCM.transform(lepPlus);
212        lepMinus = ttCM.transform(lepMinus);
213
214        // Now boost the leptons to their parent top CM frames
215        LorentzTransform topPlus, topMinus;
216        topPlus.setBetaVec(-topPlus_p4.betaVec());
217        topMinus.setBetaVec(-topMinus_p4.betaVec());
218
219        lepPlus = topPlus.transform(lepPlus);
220        lepMinus = topMinus.transform(lepMinus);
221
222        const double lepPlus_costheta_temp = lepPlus.vector3().dot(topPlus_p4.vector3()) / (lepPlus.vector3().mod() * topPlus_p4.vector3().mod());
223        const double lepMinus_costheta_temp = lepMinus.vector3().dot(topMinus_p4.vector3()) / (lepMinus.vector3().mod() * topMinus_p4.vector3().mod());
224        const double c1c2_temp = lepPlus_costheta_temp * lepMinus_costheta_temp;
225        const double cos_opening_angle_temp = lepPlus.vector3().dot(lepMinus.vector3()) / (lepPlus.vector3().mod() * lepMinus.vector3().mod());
226
227        // Fill parton-level histos
228        fillWithUFOF( _h_dphi, dphi_temp );
229        fillWithUFOF( _h_cos_opening_angle, cos_opening_angle_temp );
230        fillWithUFOF( _h_c1c2, c1c2_temp );
231        fillWithUFOF( _h_lep_costheta, lepPlus_costheta_temp );
232        fillWithUFOF( _h_lep_costheta, lepMinus_costheta_temp );
233        fillWithUFOF( _h_lep_costheta_CPV, lepPlus_costheta_temp );
234        fillWithUFOF( _h_lep_costheta_CPV, -lepMinus_costheta_temp );
235
236        // Now fill the same variables in the 2D and profile histos vs ttbar invariant mass, pT, and absolute rapidity
237        for (int i_var = 0; i_var < 3; ++i_var) {
238          double var;
239          if ( i_var == 0 ) {
240            var = tt_mass_temp;
241          } else if ( i_var == 1 ) {
242            var = tt_pT_temp;
243          } else {
244            var = tt_absrapidity_temp;
245          }
246
247          fillWithUFOF( _h_dphi_var[i_var], dphi_temp, var );
248          fillWithUFOF( _h_cos_opening_angle_var[i_var], cos_opening_angle_temp, var );
249          fillWithUFOF( _h_c1c2_var[i_var], c1c2_temp, var );
250          fillWithUFOF( _h_lep_costheta_var[i_var], lepPlus_costheta_temp, var );
251          fillWithUFOF( _h_lep_costheta_var[i_var], lepMinus_costheta_temp, var );
252          fillWithUFOF( _h_lep_costheta_CPV_var[i_var], lepPlus_costheta_temp, var );
253          fillWithUFOF( _h_lep_costheta_CPV_var[i_var], -lepMinus_costheta_temp, var );
254
255          fillWithUFOF( _h_dphi_profile[i_var], dphi_temp, var, (_h_dphi->xMax() + _h_dphi->xMin())/2. );
256          fillWithUFOF( _h_cos_opening_angle_profile[i_var], cos_opening_angle_temp, var, (_h_cos_opening_angle->xMax() + _h_cos_opening_angle->xMin())/2. );
257          fillWithUFOF( _h_c1c2_profile[i_var], c1c2_temp, var, (_h_c1c2->xMax() + _h_c1c2->xMin())/2. );
258          fillWithUFOF( _h_lep_costheta_profile[i_var], lepPlus_costheta_temp, var, (_h_lep_costheta->xMax() + _h_lep_costheta->xMin())/2. );
259          fillWithUFOF( _h_lep_costheta_profile[i_var], lepMinus_costheta_temp, var, (_h_lep_costheta->xMax() + _h_lep_costheta->xMin())/2. );
260          fillWithUFOF( _h_lep_costheta_CPV_profile[i_var], lepPlus_costheta_temp, var, (_h_lep_costheta_CPV->xMax() + _h_lep_costheta_CPV->xMin())/2. );
261          fillWithUFOF( _h_lep_costheta_CPV_profile[i_var], -lepMinus_costheta_temp, var, (_h_lep_costheta_CPV->xMax() + _h_lep_costheta_CPV->xMin())/2. );
262
263        }
264
265      }
266
267    }
268
269
270    /// Normalise histograms to unit area
271    void finalize() {
272
273      normalize(_h_dphidressedleptons);
274
275      normalize(_h_dphi);
276      normalize(_h_cos_opening_angle);
277      normalize(_h_c1c2);
278      normalize(_h_lep_costheta);
279      normalize(_h_lep_costheta_CPV);
280
281      for (int i_var = 0; i_var < 3; ++i_var) {
282        normalize(_h_dphi_var[i_var]);
283        normalize(_h_cos_opening_angle_var[i_var]);
284        normalize(_h_c1c2_var[i_var]);
285        normalize(_h_lep_costheta_var[i_var]);
286        normalize(_h_lep_costheta_CPV_var[i_var]);
287      }
288
289    }
290
291
292  private:
293
294    Histo1DPtr _h_dphidressedleptons, _h_dphi, _h_lep_costheta, _h_lep_costheta_CPV, _h_c1c2, _h_cos_opening_angle;
295    Histo2DPtr _h_dphi_var[3], _h_lep_costheta_var[3], _h_lep_costheta_CPV_var[3], _h_c1c2_var[3], _h_cos_opening_angle_var[3];
296    Profile1DPtr _h_dphi_profile[3], _h_lep_costheta_profile[3], _h_lep_costheta_CPV_profile[3], _h_c1c2_profile[3], _h_cos_opening_angle_profile[3];
297
298    const vector<double> _bins_tt_mass = {300., 430., 530., 1200.};
299    const vector<double> _bins_tt_pT = {0., 41., 92., 300.};
300    const vector<double> _bins_tt_absrapidity = {0., 0.34, 0.75, 1.5};
301    const vector<double> _bins_dphi = {0., 5.*M_PI/60., 10.*M_PI/60., 15.*M_PI/60., 20.*M_PI/60., 25.*M_PI/60., 30.*M_PI/60., 35.*M_PI/60., 40.*M_PI/60., 45.*M_PI/60., 50.*M_PI/60., 55.*M_PI/60., M_PI};
302    const vector<double> _bins_lep_costheta = {-1., -2./3., -1./3., 0., 1./3., 2./3., 1.};
303    const vector<double> _bins_lep_costheta_CPV = {-1., -2./3., -1./3., 0., 1./3., 2./3., 1.};
304    const vector<double> _bins_c1c2 = {-1., -0.4, -10./60., 0., 10./60., 0.4, 1.};
305    const vector<double> _bins_cos_opening_angle = {-1., -2./3., -1./3., 0., 1./3., 2./3., 1.};
306
307    void fillWithUFOF(Histo1DPtr h, double x) {
308      h->fill(std::max(std::min(x, h->xMax()-1e-9),h->xMin()+1e-9));
309    }
310
311    void fillWithUFOF(Histo2DPtr h, double x, double y) {
312      h->fill(std::max(std::min(x, h->xMax()-1e-9),h->xMin()+1e-9), std::max(std::min(y, h->yMax()-1e-9),h->yMin()+1e-9));
313    }
314
315    void fillWithUFOF(Profile1DPtr h, double x, double y, double c) {
316      h->fill(std::max(std::min(y, h->xMax()-1e-9),h->xMin()+1e-9), float(x > c) - float(x < c));
317    }
318
319
320  };
321
322
323  RIVET_DECLARE_PLUGIN(CMS_2016_I1413748);
324
325
326}