Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Math/Constants.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Projections/NonHadronicFinalState.hh"
#include "Rivet/Projections/VetoedFinalState.hh"
#include "Rivet/Tools/Random.hh"

namespace Rivet {


  class ATLAS_2011_I914491 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_I914491);


    /// @name Analysis methods
    /// @{

    /// Book histograms and initialise projections before the run
    void init() {

      // get the non-hadronic final-state particles
      double etaMax = 2.5;
      const NonHadronicFinalState nhfs((Cuts::etaIn(-etaMax,etaMax) && Cuts::pT >= 13.*GeV));
      declare(nhfs,"NHFS");
      // select the charged ones
      const ChargedFinalState cfs(nhfs);
      declare(cfs,"CFS");
      // and then veto electrons, and taus to be safe
      VetoedFinalState vfs(cfs);
      vfs.addVetoPairId(PID::ELECTRON);

      declare(vfs,"VFS");

      /// Book histograms
      book(_count_trigger   ,"count_trigger"  , 1, 0., 1.);
      book(_count_event     ,"count_selection", 1, 0., 1.);
      book(_count_quality   ,"count_quality"  , 1, 0., 1.);
      book(_count_beta      ,"count_beta"     , 1, 0., 1.);
      book(_count_90  ,"count_90" , 1, 0., 1.);
      book(_count_110 ,"count_110", 1, 0., 1.);
      book(_count_120 ,"count_120", 1, 0., 1.);
      book(_count_130 ,"count_130", 1, 0., 1.);

      book(_hist_beta ,"beta",1000, 0.,   2.);
      book(_hist_time ,"time",1000, -50,  50.);
      book(_hist_mass ,"mass",  60, 5., 305.);
    }


    double rndGauss(double sigma, double mean) {
      double r = sqrt(-2.0*log(rand01()));
      double phi = rand01()*2.0*pi;
      return mean+sigma*r*sin(phi);
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // smearing parameters
      // time measurement (in ns)
//       const double tsmear=5.*0.7;
      const double tsmear=0.7;
      // sagita error
      const double csag  =1.1E-4;
      // multiple scattering
      const double cms   =2.0E-2;
      // muon chamber radius (in metres)
      const double radius = 10.e3;
      // convert to ns
      const double tr = radius/c_light;
      // get the charged final-state particles
      Particles charged =
        apply<VetoedFinalState>(event,"VFS").particles();
      // need at least two candidates
      if(charged.size()<2) vetoEvent;
      // number passing trigger
      _count_trigger->fill(0.5);
      // Z mass veto
      for ( const Particle & mu1 : charged ) {
        for ( const Particle & mu2 : charged ) {
          double mass = (mu1.momentum()+mu2.momentum()).mass();
          double diff = abs(mass-91.18);
          if(diff<10.) vetoEvent;
        }
      }
      // number passing first event selection
      _count_event->fill(0.5);
      // now find the candidates
      // loop over the particles and find muons and heavy charged particles
      map<double,Particle> muonCandidates;
      for (const Particle& mu : charged) {
        // calculate the smeared momentum
        double pT     = mu.pT2();
        double pmag   = sqrt(pT+sqr(mu.pz()));
        double deltap = sqrt( sqr(csag*sqr(pmag)) + sqr(cms*mu.E()/GeV));
        double psmear = rndGauss(deltap,pmag);
        // keep particles with pT>40
        if(psmear/pmag*mu.perp()<40.*GeV||
           psmear/pmag*mu.perp()>1000.*GeV) continue;
        muonCandidates.insert(make_pair(psmear,mu));
      }
      // require two candidates
      if(muonCandidates.size()<2) vetoEvent;
      // number passing "quality" cut
      _count_quality->fill(0.5);
      // now do the time of flight
      bool filled = false;
      for(map<double,Particle>::const_iterator it=muonCandidates.begin();
          it!=muonCandidates.end();++it) {
        // true magnitude and pT of momentum
        double pT     = it->second.pT2();
        double pmag   = sqrt(pT+sqr(it->second.pz()));
        pT = sqrt(pT);
        // true time difference in ns
        double deltaT  =tr *(it->second.E()-pmag)/pT;
        // smear it
        deltaT = rndGauss(tsmear,deltaT);
        // beta
        double beta = 1./(1.+deltaT/tr*pT/pmag);
        _hist_beta->fill(beta);
        _hist_time->fill(deltaT);
        // beta cut
        if(beta<0.95) continue;
        // mass
        double mass = 2.*pT*it->first*deltaT/tr*(1.+0.5*deltaT/tr*pT/pmag);
        if(mass<0.) continue;
        mass = sqrt(mass);
        filled = true;
        _hist_mass->fill(mass);
        if(mass>90. ) {
          _count_90 ->fill(0.5);
          if(mass>110.) {
            _count_110->fill(0.5);
            if(mass>120.) {
              _count_120->fill(0.5);
              if(mass>130.) {
                _count_130->fill(0.5);
              }
            }
          }
        }
      }
      if(!filled) vetoEvent;
      // number passing beta cut
      _count_beta->fill(0.5);
    }

    /// @}

    void finalize() {
      double fact = crossSection()/picobarn/sumOfWeights()*37;
      MSG_WARNING("testing " << crossSection()/picobarn << " " << sumOfWeights() << " " << fact);
      scale(_hist_beta,fact);
      scale(_hist_time,fact);
      scale(_hist_mass,fact);
      scale( _count_trigger, fact);
      scale( _count_event, fact);
      scale( _count_quality, fact);
      scale( _count_beta, fact);
      scale( _count_90, fact);
      scale( _count_110, fact);
      scale( _count_120, fact);
      scale( _count_130, fact);
    }

  private:

    /// @name Histograms
    /// @{
    Histo1DPtr _hist_beta;
    Histo1DPtr _hist_time;
    Histo1DPtr _hist_mass;
    Histo1DPtr _count_trigger;
    Histo1DPtr _count_event;
    Histo1DPtr _count_quality;
    Histo1DPtr _count_beta;
    Histo1DPtr _count_90;
    Histo1DPtr _count_110;
    Histo1DPtr _count_120;
    Histo1DPtr _count_130;
    /// @}

  };


  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_I914491, ATLAS_2011_S9108483);

}