Ticket #74: MYANALYSIS.cc

#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/AnalysisLoader.hh"

namespace Rivet {
  

  class MYANALYSIS : public Analysis {
  public:

    MYANALYSIS() : Analysis("MYANALYSIS") {}

    void init() {
      IdentifiedFinalState piplus;
      piplus.acceptId(PIPLUS);
      addProjection(piplus, "piplus");
      
      IdentifiedFinalState piminus;
      piminus.acceptId(PIMINUS);
      addProjection(piminus, "piminus");
      
      IdentifiedFinalState neutrinos;
      neutrinos.acceptNeutrinos();
      addProjection(neutrinos, "neutrinos");
    }
    
    void analyze(const Event& event) {
      ParticleVector piplusv = applyProjection<FinalState>(event, "piplus").particles();
      ParticleVector piminusv = applyProjection<FinalState>(event, "piminus").particles();
      ParticleVector neutrinosv = applyProjection<FinalState>(event, "neutrinos").particles();
      
      if (piplusv.size()!=1) {
        std::cout<<"No piplus"<<std::endl;
        vetoEvent;
      }
      if (piminusv.size()!=1) {
        std::cout<<"No piminus"<<std::endl;
        vetoEvent;
      }
      if (neutrinosv.size()!=2) {
        std::cout<<"No neutrinos"<<std::endl;
        vetoEvent;
      }
      
      FourMomentum restframe(piplusv[0].momentum()+piminusv[0].momentum()+neutrinosv[0].momentum()+neutrinosv[1].momentum());
      
      LorentzTransform cms_boost(-restframe.boostVector());
      FourMomentum piplus(cms_boost.transform(piplusv[0].momentum()));
      FourMomentum piminus(cms_boost.transform(piminusv[0].momentum()));
      
      double zplus = 2.0*piplus.E()/120.0;
      double zminus = 2.0*piminus.E()/120.0;
      
      if (zplus>0.5 && zminus>0.5) {
        _sigma1+=event.weight();
      }
      else if (zplus<0.5 && zminus<0.5) {
        _sigma2+=event.weight();
      }
      else if (zplus>0.5 && zminus<0.5) {
        _sigma3+=event.weight();
      }
      else if (zplus<0.5 && zminus>0.5) {
        _sigma4+=event.weight();
      }
    }
    
    void finalize() {
      // No histos, so nothing to do!
      std::cout<<"_sigma1="<<_sigma1*crossSectionPerEvent()<<std::endl;
      std::cout<<"_sigma2="<<_sigma2*crossSectionPerEvent()<<std::endl;
      std::cout<<"_sigma3="<<_sigma3*crossSectionPerEvent()<<std::endl;
      std::cout<<"_sigma4="<<_sigma4*crossSectionPerEvent()<<std::endl;
      
      std::cout<<"A_fastslow="<<(_sigma1+_sigma2-_sigma3-_sigma4)/(_sigma1+_sigma2+_sigma3+_sigma4)<<std::endl;
    }
    
  private:
    
    double _sigma1, _sigma2, _sigma3, _sigma4;
  };
  
  

  // This global object acts as a hook for the plugin system
  AnalysisBuilder<MYANALYSIS> plugin_MYANALYSIS;
  
}