Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/ChargedFinalState.hh"

namespace Rivet {

  /// @brief Study of ordered hadron chains at 7 TeV
  class ATLAS_2017_I1624693 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1624693);

    /// @name Analysis methods
    /// @{
    struct usedX {

      int locMin;
      int locMax;
      std::vector<std::pair<int,float> > chains;

      // Constructor
      usedX(int min, int max, int ic, float mass) {
        locMin=min;
        locMax=max;
        chains.clear();
        chains.push_back(std::pair<int,float>(ic,mass));
      }

      // Constructor
      usedX(int min, int max) {
        locMin=min;
        locMax=max;
        chains.clear();
      }

      void add(int jc, float mass) {

        if (chains.size()) {
          std::vector<std::pair<int,float> >::iterator it=chains.begin();
          while ( it!=chains.end() && mass>(*it).second )  ++it;
          chains.insert(it,std::pair<int,float>(jc,mass));
        }
        else {
          chains.push_back(std::pair<int,float>(jc,mass));
        }
      }
    };


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

      /// @todo Initialise and register projections here
      ChargedFinalState cfs((Cuts::etaIn(-2.5, 2.5) && Cuts::pT >=  0.1*GeV));
      declare(cfs,"CFS");

      // pion mass;
      pim = 0.1396;

      /// @todo Book histograms here, e.g.:
      book(_DeltaQ , 1, 1, 1);
      book(_Delta3h, 2, 1, 1);
      book(_dalitz , 3, 1, 1);

      // auxiliary
      book(_h_nch, "_nch", 200, -0.5, 199.5);

    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      //const double weight = event.weight();
      bool match =false;

      /// @todo Do the event by event analysis here
      const ChargedFinalState& had = apply<ChargedFinalState>(event, "CFS");
      Particles hs=had.particles();
      int nch = hs.size();

      if (nch < 3)  return;

      _h_nch->fill(1.*nch,1.);

      for (unsigned int i=0; i < hs.size() - 1; ++i) {
        for (unsigned int j=i+1; j < hs.size(); ++j) {
          double q12 = qq(hs[i],hs[j],match);
          if (match) _DeltaQ->fill(q12,-1.);
          else       _DeltaQ->fill(q12,1.);
        }
      }

      // chain selection

      std::vector<float> wchain;
      std::vector< std::vector<unsigned int> > rchains;
      std::vector< std::vector<float> > mchains;
      wchain.clear();  rchains.clear();  mchains.clear();
      for (unsigned int ip1 = 0; ip1< hs.size(); ++ip1 ) {
        wchain.push_back(1.);
        std::vector<unsigned int> cc(1,ip1);
        std::vector<float> mc;

        double qlmin=10000.; int ilmin=-1;
        for (unsigned ip2 = 0; ip2 < hs.size(); ++ip2) {
          if (ip2==ip1) continue;
          double ql = qq(hs[ip1],hs[ip2],match);
          if (!match) continue;    // looking for closest like-sign match
          if (ql <qlmin) { qlmin=ql; ilmin=ip2;}
        }
        if (ilmin<0) {
          wchain.back()=0.;
          mc.push_back(-1.);
        }
        else {     // search for unlike-sign match
          cc.push_back(ilmin);
          mc.push_back(qlmin);
          if (int(ip1)>ilmin && rchains[ilmin][1]==ip1) {
            // std::cout <<"exclusive match:"<< std::endl;
            wchain.back()=0.5; wchain[ilmin]=0.5;
          }

          double m3min=10000.; int ixmin=-1;
          for (unsigned ip2 = 0; ip2< hs.size(); ++ip2) {
            if (ip2==ip1 || int(ip2)==ilmin )  continue;
            double qx = qq(hs[ip1],hs[ip2],match);
            if (match)  continue;
            double qxl = qq(hs[ip2],hs[ilmin],match);
            double m3 = sqrt(9*pim*pim+qxl*qxl+qlmin*qlmin+qx*qx);
            if (m3 <m3min) { m3min=m3; ixmin=ip2;}
          }

          if (ixmin<0) {
            wchain.back()=0.;
            mc.push_back(-1.);
          }
          else {
            cc.push_back(ixmin);
            mc.push_back(m3min);
          }
        }
        rchains.push_back(cc);
        mchains.push_back(mc);
      }

      // cleanup: association rate for like-sign pairs should not exceed 2
      std::vector<float> assoc(hs.size(),0.);       // cache for association rate
      std::vector<bool> accept(rchains.size(), false);
      // loop over chains and accept lowest masses while watching the association rate
      int inext = 0;
      while ( inext>-1 ) {
        inext = -1; float cMin = 100000.;
        // find non-accepted chain with lowest Q_ls; dissolve chains if association count over 2
        for (unsigned int ic=0; ic < rchains.size(); ++ic) {
          if (rchains[ic].size() < 2)  continue;
          if (accept[ic])  continue;
          if (mchains[ic][0] < cMin) { cMin = mchains[ic][0]; inext=ic; }
        }
        if (inext>-1 ) {
          unsigned int cloc0 = rchains[inext][0];
          unsigned int cloc1 = rchains[inext][1];
          if ( (assoc[cloc0] + 1. <= 2.)  &&  (assoc[cloc1] + 1. <= 2.) ) { // chain can be accepted
            accept[inext]=true;
            assoc[cloc0]+=1.;
            assoc[cloc1]+=1.;
            if (wchain[inext]==0.5) {  // accept the identical chain, too
              for (unsigned int ic=0; ic<hs.size(); ++ic) {
                if (rchains[ic][0] == cloc1 && rchains[ic][1] == cloc0) {
                  accept[ic]=true;
                  break;
                }
              }
            }
          }
          else if ( assoc[cloc0]>1 ) { // association count filled up, discard chain
            accept[inext]=true;
            wchain[inext]=0.;
          }
          else {   // dissolve chain and find new association
            unsigned int i1 = rchains[inext][0];
            float mMn = 1000000.;
            int ipn = -1;
            for (unsigned int i2=0; i2<hs.size(); ++i2) {
              if (i1 == i2)  continue;
              double m = qq(hs[i1],hs[i2],match);
              if (!match)  continue;
              if (assoc[i2] > 1.)  continue;
              if (m > 0. && m <mMn ) { mMn = m; ipn = i2;}
            }
            if (ipn >= 0) {
              rchains[inext][1]=ipn;  mchains[inext][0]=mMn;
              // resolve chain weight : by default, it is 1.
              wchain[inext]=1.;
              // check exclusivity of pairing
              for (unsigned int ico=0; ico<hs.size(); ++ico) {
                if (int(rchains[ico][0]) == ipn && rchains[ico][1] == i1) {   // scale the contribution from both chains
                  wchain[ico]=0.5;
                  wchain[inext]=0.5;
                }
              }
              // add 3.member
              // continue with arbitrary match
              int ipnn=-1; float mMnn = 10000.;
              mMn = 1000000.;
              for (unsigned int ij=0; ij < hs.size(); ++ij) {
                rchains[inext].resize(2);
                float q02 = qq(hs[i1],hs[ij],match);
                if (match>0.)  continue;
                float q12 = qq(hs[ipn],hs[ij],match);
                double m3 = sqrt(9*pim*pim+q02*q02+mMn*mMn+q12*q12);
                if (m3>0. && m3 <mMnn ) { mMnn = m3; ipnn = ij; }
              }
              if (ipnn>=0) { rchains[inext].push_back(ipnn); rchains[inext][2]=ipnn; mchains[inext][1]=mMnn; }
              else {accept[inext]=true; wchain[inext]=0.;}
            }
            else { // chain not recovered
              wchain[inext]=0.;
              accept[inext]=true;
            }
          }
        }
      }  // end loop over chains

      // cleanup: association rate for unlike-sign pairs
      // third member verification
      std::vector<bool> accept3(rchains.size(),false);
      // watch unlike-sign combinations used
      std::vector<usedX> used;
      // loop over chains and accept lowest masses while watching the association rate
      inext = 0;
      while ( inext>-1 ) {
        inext = -1; float cMin = 100000.;
        // find non-accepted chain with lowest mass; dissolve chains if association count over 3
        for (unsigned int ic=0; ic < rchains.size(); ++ic) {
          if (rchains[ic].size() < 3 || !wchain[ic] || !accept[ic])  continue;
          if (accept3[ic])  continue;
          if (mchains[ic][1]<cMin) { cMin = mchains[ic][1]; inext=ic; }
        }
        // check association counts
        if (inext>-1 ) {
          unsigned int cloc0 = rchains[inext][0];
          unsigned int cloc1 = rchains[inext][1];
          unsigned int cloc2 = rchains[inext][2];

          // map use of unlike sign pairs
          int iu0 = -1; float w0=0.;
          for (unsigned int iu=0; iu<used.size(); ++iu) {
            if (fmin(cloc0,cloc2)==used[iu].locMin && fmax(cloc0,cloc2)==used[iu].locMax ) {
              iu0=iu;
              if (used[iu].chains.size() > 0)
                for (unsigned int iw=0; iw<used[iu].chains.size(); ++iw)  w0+=wchain[used[iu].chains[iw].first];
              //used[iu].add(i1,mch[1]);
              break;
            }
          }
          if (iu0<0) { used.push_back(usedX(fmin(cloc0,cloc2),fmax(cloc0,cloc2)));iu0=used.size()-1; }
          int iu1 = -1; float w1=0.;
          for (unsigned int iu=0; iu<used.size(); ++iu) {
            if (fmin(cloc1,cloc2)==used[iu].locMin && fmax(cloc1,cloc2)==used[iu].locMax) {
              iu1=iu;
              if (used[iu].chains.size()>0)
                for (unsigned int iw=0; iw<used[iu].chains.size(); iw++) w1 += wchain[used[iu].chains[iw].first];
              //used[iu].add(inext,mch[1]);
              break;
            }
          }
          if (iu1<0) { used.push_back(usedX(fmin(cloc1,cloc2),fmax(cloc1,cloc2))); iu1=used.size()-1; }

          if ( assoc[cloc2] < 3. && w0 < 2. && w1 < 2.) {
            accept3[inext] = true;
            assoc[cloc2] += 1.;
            used[iu0].add(inext, mchains[inext][1]);
            used[iu1].add(inext, mchains[inext][1]);
            if (wchain[inext]==0.5) {  // accept the identical chain, too
              for (unsigned int ic=0; ic< rchains.size(); ++ic) {
                if (rchains[ic][0]==cloc1 && rchains[ic][1] == cloc0) {
                  accept3[ic]=true;
                  used[iu0].add(ic, mchains[ic][1]);
                  used[iu1].add(ic, mchains[ic][1]);
                  break;
                }
              }
            }
          }
          else { // find new association
            int i1 = rchains[inext][0];
            int i2 = rchains[inext][1];
            float mMn = 1000000.;
            int ipn=-1; int iploc=-1;
            rchains[inext].pop_back();
            for (unsigned int i3 = 0; i3 < hs.size(); ++i3) {
              double q02 = qq(hs[i1],hs[i3],match);
              if (match > 0.)  continue;
              if (assoc[i3] > 3-wchain[inext])  continue;
              // check pair association
              w0=0.; w1=0.;
              for (unsigned int iu=0; iu<used.size(); ++iu) {
                if (fmin(cloc0,i3)==used[iu].locMin && fmax(cloc0,i3)==used[iu].locMax ) {
                  if (used[iu].chains.size() > 0)
                    for (unsigned int iw=0; iw<used[iu].chains.size(); ++iw)  w0 += wchain[used[iu].chains[iw].first];
                }
                if (fmin(cloc1,i3)==used[iu].locMin && fmax(cloc1,i3)==used[iu].locMax ) {
                  if (used[iu].chains.size()>0)
                    for (unsigned int iw=0; iw<used[iu].chains.size(); ++iw)  w1 += wchain[used[iu].chains[iw].first];
                }
              }
              if (w0+wchain[inext]>2. || w1+wchain[inext]>2.) continue;

              float q12 = qq(hs[i2],hs[i3],match);
              float q01 = qq(hs[i1],hs[i2],match);
              float m = sqrt(9*pim*pim+q02*q02+q01*q01+q12*q12);
              if (m>0. && m <mMn ) { mMn = m; ipn = i3; iploc = i3; }
            }
            if (ipn>=0) {
              rchains[inext].push_back(ipn); rchains[inext][2]=iploc; mchains[inext][1]=mMn;
            }
            else { // chain not recovered
              wchain[inext]=0.;
            }
          }
        }
      }  // end loop over chains
      // end 3rd member optimization

      for (unsigned int ip=0; ip < wchain.size(); ++ip) {
        if (!wchain[ip])  continue;
        if (rchains[ip].size() < 3)  continue;
        float m3min = mchains[ip][1];
        if (m3min > 0.59)  continue;
        // dalitz plot
        std::pair<float,float> dd = dalitz3(hs[rchains[ip][0]], hs[rchains[ip][1]], hs[rchains[ip][2]]);
        _dalitz->fill(dd.first,dd.second,1.*wchain[ip]);
        // Delta(Q) spectra
        float qlmin = mchains[ip][0];
        float qxmin = qq(hs[rchains[ip][0]], hs[rchains[ip][2]], match);
        float xlmin = qq(hs[rchains[ip][1]], hs[rchains[ip][2]], match);
        _Delta3h->fill(qxmin, 0.5*wchain[ip]);
        _Delta3h->fill(xlmin, 0.5*wchain[ip]);
        _Delta3h->fill(qlmin, -1.*wchain[ip]);
      }
    }

    /// Normalise histograms etc., after the run
    void finalize() {


      // normalize by the number of charged particles
      // counter automatic division by bin size
      double norm = 0.01 / (_h_nch->xMean()*_h_nch->numEntries());
      _dalitz->scaleW(norm);
      _DeltaQ->scaleW(norm);
      _Delta3h->scaleW(norm);

    }

    /// @}
    double qq(const Particle& gp1, const Particle& gp2, bool& match) {
      match = gp1.charge() * gp2.charge() > 0;
      FourMomentum p1, p2;
      p1.setPM(gp1.px(), gp1.py(), gp1.pz(), pim);
      p2.setPM(gp2.px(), gp2.py(), gp2.pz(), pim);
      return sqrt(fmax(0., (p1 + p2).mass2() - 4*pim*pim));
    }

    std::pair<float,float> dalitz3(const Particle& gp1, const Particle& gp2, const Particle& gp3) const {

      float p1= gp1.pt();
      float p2= gp2.pt();
      float p3= gp3.pt();
      float th1 = gp1.theta();
      float th2 = gp2.theta();
      float th3 = gp3.theta();
      float ph1 = gp1.phi();
      float ph2 = gp2.phi();
      float ph3 = gp3.phi();
      float e1 = sqrt(p1*p1+pim*pim);
      float e2 = sqrt(p2*p2+pim*pim);
      float e3 = sqrt(p3*p3+pim*pim);

      float p1x = p1*cos(ph1)*sin(th1);
      float p1y = p1*sin(ph1)*sin(th1);
      float p1z = p1*cos(th1);

      float p2x = p2*cos(ph2)*sin(th2);
      float p2y = p2*sin(ph2)*sin(th2);
      float p2z = p2*cos(th2);

      float p3x = p3*cos(ph3)*sin(th3);
      float p3y = p3*sin(ph3)*sin(th3);
      float p3z = p3*cos(th3);

      float px = p1x+p2x+p3x;
      float py = p1y+p2y+p3y;
      float pz = p1z+p2z+p3z;
      float ap = sqrt(px*px+py*py+pz*pz);
      float e=e1+e2+e3;

      float beta = ap/e;
      float gamma = 1./sqrt(1-beta*beta);

      float p1l = (p1x*px+p1y*py+p1z*pz)/ap;
      float p2l = (p2x*px+p2y*py+p2z*pz)/ap;
      float p3l = (p3x*px+p3y*py+p3z*pz)/ap;

      float e1_boost = gamma*e1-gamma*beta*p1l;
      float e2_boost = gamma*e2-gamma*beta*p2l;
      float e3_boost = gamma*e3-gamma*beta*p3l;

      float Q = sqrt(e*e-ap*ap)-3*pim;

      return std::pair<float,float>(sqrt(3.)*(e1_boost-e2_boost)/Q , 3*(e3_boost-pim)/Q-1.);
    }

  private:

    // Data members like post-cuts event weight counters go here
    float pim;

  private:

    /// @name Histograms
    Histo1DPtr  _DeltaQ;
    Histo1DPtr _Delta3h;
    Histo1DPtr   _h_nch;
    Histo2DPtr _dalitz;
    /// @}
  };

  RIVET_DECLARE_PLUGIN(ATLAS_2017_I1624693);
}