1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/Thrust.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief p pbar correlations
 10  class DELPHI_2000_I531568 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_2000_I531568);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Initialise and register projections
 24      const ChargedFinalState cfs;
 25      declare(cfs        ,"CFS"   );
 26      declare(Thrust(cfs),"Thrust");
 27      // book histos
 28      book(_h_pMp,"_n_pMp",8,0.,1.);
 29      book(_h_sum,"_n_sum",8,0.,1.);
 30    }
 31
 32    void findPP(unsigned int mode, const Vector3 & axis, const Particles & part,
 33		unsigned int & pp, double & dy) {
 34      pp=0;
 35      dy=1e30;
 36      map<double,Particle> rapOrdered;
 37      for(const Particle & p : part) {
 38      	const double mom = dot(axis, p.momentum().p3());
 39      	const double energy = p.E();
 40      	const double rap = 0.5 * std::log((energy + mom) / (energy - mom));
 41      	if(mode==0) {
 42      	  if(rap>0.)  rapOrdered[rap] = p;
 43      	}
 44      	else {
 45      	  if(rap<=0.) rapOrdered[rap] = p;
 46      	}
 47      }
 48      map<double,Particle>::const_iterator ii[2]={rapOrdered.end(),rapOrdered.end()};
 49      // map<double,Particle>::const_iterator it,i,im;
 50      for(map<double,Particle>::const_iterator it=rapOrdered.begin();it!=rapOrdered.end();++it) {
 51      	if(it->second.abspid()==PID::PROTON) {
 52	  if(ii[0]!=rapOrdered.end())
 53	    ii[1]=it;
 54	  else
 55	    ii[0]=it;
 56	}
 57      }
 58      // number of particles between the proton and antiproton
 59      int rank = std::distance(ii[0],ii[1]);
 60      if(rank>4) return;
 61      // protom/antiproton next to each other, distance to nearest meson near them
 62      if(rank==1) {
 63	map<double,Particle>::const_iterator im=ii[0];--ii[0];
 64	map<double,Particle>::const_iterator ip=ii[1];++ii[1];
 65	if(ii[0]!=rapOrdered.begin()) {
 66	  pp=1;
 67	  dy = min(2./3.*abs(im->first-ii[0]->first),dy);
 68	}
 69	if(ip!=rapOrdered.end()) {
 70	  pp=1;
 71	  dy = min(2./3.*abs(ip->first-ii[1]->first),dy);
 72	}
 73      }
 74      else {
 75	double ycent = 0.5*(ii[0]->first+ii[1]->first);
 76	double ymin=1e30;
 77	map<double,Particle>::const_iterator im=rapOrdered.end();
 78	map<double,Particle>::const_iterator it=ii[0];++it;
 79	for(;it!=ii[1];++it) {
 80	  double test = abs(ycent-it->first);
 81	  if(test<ymin) {
 82	    im=it;
 83	    ymin=test;
 84	  }
 85	}
 86	pp=2;
 87	dy = min(abs(ii[0]->first-im->first),abs(ii[1]->first-im->first));
 88      }
 89    }
 90
 91    /// Perform the per-event analysis
 92    void analyze(const Event& event) {
 93      const Thrust thrust = apply<Thrust>(event,"Thrust");
 94      Vector3 axis = thrust.thrustAxis();
 95      const ChargedFinalState cfs = apply<ChargedFinalState>(event,"CFS");
 96      unsigned int np[2]={0,0}, npbar[2]={0,0};
 97      for(const Particle & p : cfs.particles()) {
 98        const double mom = dot(axis, p.momentum().p3());
 99        const double energy = p.E();
100        const double rap = 0.5 * std::log((energy + mom) / (energy - mom));
101	if(p.abspid()==PID::PROTON) {
102	  unsigned int irap = rap>0 ? 0 : 1;
103	  if(p.pid()>0) ++np   [irap];
104	  else          ++npbar[irap];
105	}
106      }
107      if(np[0]==1 && npbar[0]==1) {
108	unsigned int pp = 0;
109	double dy(1e30);
110	findPP(0,axis,cfs.particles(),pp,dy);
111	if(pp==1) {
112	  _h_sum->fill(dy);
113	}
114	else if(pp==2) {
115	  _h_sum->fill(dy);
116	  _h_pMp->fill(dy);
117	}
118      }
119      if(np[1]==1 && npbar[1]==1) {
120	unsigned int pp = 0;
121	double dy(1e30);
122	findPP(1,axis,cfs.particles(),pp,dy);
123	if(pp==1) {
124	  _h_sum->fill(dy);
125	}
126	else if(pp==2) {
127	  _h_sum->fill(dy);
128	  _h_pMp->fill(dy);
129	}
130      }
131    }
132
133
134    /// Normalise histograms etc., after the run
135    void finalize() {
136      Estimate1DPtr h_r;
137      book(h_r,1,1,1);
138      divide(_h_pMp,_h_sum,h_r);
139    }
140
141    /// @}
142
143
144    /// @name Histograms
145    /// @{
146    Histo1DPtr _h_sum,_h_pMp;
147    /// @}
148
149
150  };
151
152
153  RIVET_DECLARE_PLUGIN(DELPHI_2000_I531568);
154
155}