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#include "Rivet/Tools/Random.hh"
  8
  9namespace Rivet {
 10
 11
 12  class ATLAS_2011_I914491 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_I914491);
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25      // get the non-hadronic final-state particles
 26      double etaMax = 2.5;
 27      const NonHadronicFinalState nhfs((Cuts::etaIn(-etaMax,etaMax) && Cuts::pT >= 13.*GeV));
 28      declare(nhfs,"NHFS");
 29      // select the charged ones
 30      const ChargedFinalState cfs(nhfs);
 31      declare(cfs,"CFS");
 32      // and then veto electrons, and taus to be safe
 33      VetoedFinalState vfs(cfs);
 34      vfs.addVetoPairId(PID::ELECTRON);
 35
 36      declare(vfs,"VFS");
 37
 38      /// Book histograms
 39      book(_count_trigger   ,"count_trigger"  , 1, 0., 1.);
 40      book(_count_event     ,"count_selection", 1, 0., 1.);
 41      book(_count_quality   ,"count_quality"  , 1, 0., 1.);
 42      book(_count_beta      ,"count_beta"     , 1, 0., 1.);
 43      book(_count_90  ,"count_90" , 1, 0., 1.);
 44      book(_count_110 ,"count_110", 1, 0., 1.);
 45      book(_count_120 ,"count_120", 1, 0., 1.);
 46      book(_count_130 ,"count_130", 1, 0., 1.);
 47
 48      book(_hist_beta ,"beta",1000, 0.,   2.);
 49      book(_hist_time ,"time",1000, -50,  50.);
 50      book(_hist_mass ,"mass",  60, 5., 305.);
 51    }
 52
 53
 54    double rndGauss(double sigma, double mean) {
 55      double r = sqrt(-2.0*log(rand01()));
 56      double phi = rand01()*2.0*pi;
 57      return mean+sigma*r*sin(phi);
 58    }
 59
 60    /// Perform the per-event analysis
 61    void analyze(const Event& event) {
 62      // smearing parameters
 63      // time measurement (in ns)
 64//       const double tsmear=5.*0.7;
 65      const double tsmear=0.7;
 66      // sagita error
 67      const double csag  =1.1E-4;
 68      // multiple scattering
 69      const double cms   =2.0E-2;
 70      // muon chamber radius (in metres)
 71      const double radius = 10.e3;
 72      // convert to ns
 73      const double tr = radius/c_light;
 74      // get the charged final-state particles
 75      Particles charged =
 76        apply<VetoedFinalState>(event,"VFS").particles();
 77      // need at least two candidates
 78      if(charged.size()<2) vetoEvent;
 79      // number passing trigger
 80      _count_trigger->fill(0.5);
 81      // Z mass veto
 82      for ( const Particle & mu1 : charged ) {
 83        for ( const Particle & mu2 : charged ) {
 84          double mass = (mu1.momentum()+mu2.momentum()).mass();
 85          double diff = abs(mass-91.18);
 86          if(diff<10.) vetoEvent;
 87        }
 88      }
 89      // number passing first event selection
 90      _count_event->fill(0.5);
 91      // now find the candidates
 92      // loop over the particles and find muons and heavy charged particles
 93      map<double,Particle> muonCandidates;
 94      for (const Particle& mu : charged) {
 95        // calculate the smeared momentum
 96        double pT     = mu.pT2();
 97        double pmag   = sqrt(pT+sqr(mu.pz()));
 98        double deltap = sqrt( sqr(csag*sqr(pmag)) + sqr(cms*mu.E()/GeV));
 99        double psmear = rndGauss(deltap,pmag);
100        // keep particles with pT>40
101        if(psmear/pmag*mu.perp()<40.*GeV||
102           psmear/pmag*mu.perp()>1000.*GeV) continue;
103        muonCandidates.insert(make_pair(psmear,mu));
104      }
105      // require two candidates
106      if(muonCandidates.size()<2) vetoEvent;
107      // number passing "quality" cut
108      _count_quality->fill(0.5);
109      // now do the time of flight
110      bool filled = false;
111      for(map<double,Particle>::const_iterator it=muonCandidates.begin();
112          it!=muonCandidates.end();++it) {
113        // true magnitude and pT of momentum
114        double pT     = it->second.pT2();
115        double pmag   = sqrt(pT+sqr(it->second.pz()));
116        pT = sqrt(pT);
117        // true time difference in ns
118        double deltaT  =tr *(it->second.E()-pmag)/pT;
119        // smear it
120        deltaT = rndGauss(tsmear,deltaT);
121        // beta
122        double beta = 1./(1.+deltaT/tr*pT/pmag);
123        _hist_beta->fill(beta);
124        _hist_time->fill(deltaT);
125        // beta cut
126        if(beta<0.95) continue;
127        // mass
128        double mass = 2.*pT*it->first*deltaT/tr*(1.+0.5*deltaT/tr*pT/pmag);
129        if(mass<0.) continue;
130        mass = sqrt(mass);
131        filled = true;
132        _hist_mass->fill(mass);
133        if(mass>90. ) {
134          _count_90 ->fill(0.5);
135          if(mass>110.) {
136            _count_110->fill(0.5);
137            if(mass>120.) {
138              _count_120->fill(0.5);
139              if(mass>130.) {
140                _count_130->fill(0.5);
141              }
142            }
143          }
144        }
145      }
146      if(!filled) vetoEvent;
147      // number passing beta cut
148      _count_beta->fill(0.5);
149    }
150
151    /// @}
152
153    void finalize() {
154      double fact = crossSection()/picobarn/sumOfWeights()*37;
155      MSG_WARNING("testing " << crossSection()/picobarn << " " << sumOfWeights() << " " << fact);
156      scale(_hist_beta,fact);
157      scale(_hist_time,fact);
158      scale(_hist_mass,fact);
159      scale( _count_trigger, fact);
160      scale( _count_event, fact);
161      scale( _count_quality, fact);
162      scale( _count_beta, fact);
163      scale( _count_90, fact);
164      scale( _count_110, fact);
165      scale( _count_120, fact);
166      scale( _count_130, fact);
167    }
168
169  private:
170
171    /// @name Histograms
172    /// @{
173    Histo1DPtr _hist_beta;
174    Histo1DPtr _hist_time;
175    Histo1DPtr _hist_mass;
176    Histo1DPtr _count_trigger;
177    Histo1DPtr _count_event;
178    Histo1DPtr _count_quality;
179    Histo1DPtr _count_beta;
180    Histo1DPtr _count_90;
181    Histo1DPtr _count_110;
182    Histo1DPtr _count_120;
183    Histo1DPtr _count_130;
184    /// @}
185
186  };
187
188
189  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_I914491, ATLAS_2011_S9108483);
190
191}