Ticket #74: MYANALYSIS.2.cc

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

namespace Rivet {


  class MYANALYSIS : public Analysis {
    public:

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

      void init() {
        _sigma1 = 0;
        _sigma2 = 0;
        _sigma3 = 0;
        _sigma4 = 0;
      }

      void analyze(const Event& event) {

        GenEvent::particle_const_iterator p;
        vector<vector <GenParticle> > tau_parts;
        vector<GenParticle> tmp_parts;
        //loop over all particles
        for (p = event.genEvent().particles_begin();
            p != event.genEvent().particles_end(); ++p) {
          //loop over particles at decay vertex of status = 2 taus 
          if  (fabs((*p)->pdg_id())==15 && (*p)->status()==2){
            const GenVertex* dv = (*p)->end_vertex();

            bool found = 0;
            int count=0;
            tmp_parts.clear();
            for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin() ;
                pp != dv->particles_out_const_end(); ++pp) {
              if (fabs((*pp)->pdg_id())==16) found=1;
              tmp_parts.push_back(**pp);
              count++;
            }
            //are there exactly two particles at this tau decay vertex, one of which is a tau (anti) neutrino 
            if (count !=2) found =0;
            if (found){
              //then push back these two particles as a vector
              tau_parts.push_back(tmp_parts);
            } 
          } 
        }
        //are there exactly two vectors of two particles, with a tau (anti)neutrino in each?
        if (tau_parts.size()!=2) vetoEvent;


        GenParticle *piplus, *piminus;

        float re, rpx, rpy, rpz;
        re = rpx = rpy = rpz = 0;
        bool find_tauplus,find_tauminus;
        find_tauplus = find_tauminus = 0;

        //cout <<"ntau = "<<tau_parts.size()<<": children are\t";  
        for (unsigned int i=0;i<tau_parts.size();i++){
          for (unsigned int j=0;j<tau_parts.at(i).size();j++){
            //cout <<"\t\t"<<tau_parts.at(i).at(j).pdg_id();
            if (tau_parts.at(i).at(j).pdg_id()==211){  piplus=  &GenParticle(tau_parts.at(i).at(j));find_tauplus=1;}
            if (tau_parts.at(i).at(j).pdg_id()==-211){ piminus= &GenParticle(tau_parts.at(i).at(j));find_tauminus=1;}
            re += tau_parts.at(i).at(j).momentum().e();
            rpx += tau_parts.at(i).at(j).momentum().px();
            rpy += tau_parts.at(i).at(j).momentum().py();
            rpz += tau_parts.at(i).at(j).momentum().pz();
          }
          //cout <<"\n";
        }
        //cout <<"\n";

        //failed to find pi+-, somehow
        if (!find_tauplus || !find_tauminus) vetoEvent;

        //make calculation
        FourMomentum restframe(re,rpx,rpy,rpz);
        LorentzTransform cms_boost(-restframe.boostVector());
        FourMomentum piplus_b(cms_boost.transform(piplus->momentum()));
        FourMomentum piminus_b(cms_boost.transform(piminus->momentum()));

        double zplus = 2.0*piplus_b.E()/120.0;
        double zminus = 2.0*piminus_b.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();
        }

        if (zplus <0 || zminus <0 || zplus>1. || zplus>1.){
          std::cout <<"zplus = "<<zplus<<"\t zminus = "<<zminus<<std::endl;
        }
      }

      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;

        std::cout<<"_sigma1="<<_sigma1<<std::endl;
        std::cout<<"_sigma2="<<_sigma2<<std::endl;
        std::cout<<"_sigma3="<<_sigma3<<std::endl;
        std::cout<<"_sigma4="<<_sigma4<<std::endl;

      }

    private:

      double _sigma1, _sigma2, _sigma3, _sigma4;
  };



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


}