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}