1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/UnstableParticles.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief e+e- > rho+ rho-
 11  class BABAR_2008_I789011 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2008_I789011);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23      // Initialise and register projections
 24      declare(Beam(), "Beams");
 25      declare(UnstableParticles(Cuts::abspid==213), "UFS");
 26      declare(FinalState(), "FS");
 27      // histos
 28      for(unsigned int ix=0;ix<5;++ix)
 29	book(_h[ix],1,1,1+ix);
 30    }
 31
 32    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 33      for( const Particle &child : p.children()) {
 34	if(child.children().empty()) {
 35	  nRes[child.pid()]-=1;
 36	  --ncount;
 37	}
 38	else
 39	  findChildren(child,nRes,ncount);
 40      }
 41    }
 42
 43    /// Perform the per-event analysis
 44    void analyze(const Event& event) {
 45      // get the axis, direction of incoming electron
 46      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 47      Vector3 axis;
 48      if(beams.first.pid()>0)
 49	axis = beams.first .momentum().p3().unit();
 50      else
 51	axis = beams.second.momentum().p3().unit();
 52      // types of final state particles
 53      const FinalState& fs = apply<FinalState>(event, "FS");
 54      map<long,int> nCount;
 55      int ntotal(0);
 56      for (const Particle& p :  fs.particles()) {
 57	nCount[p.pid()] += 1;
 58	++ntotal;
 59      }
 60      // loop over rho mesons
 61      const UnstableParticles & ufs = apply<UnstableParticles>(event, "UFS");
 62      Particle rhoP,rhoM;
 63      bool matched(false);
 64      for (const Particle& p :  ufs.particles()) {
 65	if(p.children().empty()) continue;
 66	map<long,int> nRes=nCount;
 67	int ncount = ntotal;
 68	findChildren(p,nRes,ncount);
 69	matched=false;
 70	// check for antiparticle
 71	for (const Particle& p2 :  ufs.particles(Cuts::pid==-p.pid())) {
 72	  if(p2.children().empty()) continue;
 73	  map<long,int> nRes2=nRes;
 74	  int ncount2 = ncount;
 75	  findChildren(p2,nRes2,ncount2);
 76	  if(ncount2==0) {
 77	    matched = true;
 78	    for(auto const & val : nRes2) {
 79	      if(val.second!=0) {
 80		matched = false;
 81		break;
 82	      }
 83	    }
 84	    // found rho+/-
 85	    if(matched) {
 86	      if(p.pid()>0) {
 87		rhoP = p;
 88		rhoM = p2;
 89	      }
 90	      else {
 91		rhoP = p2;
 92		rhoM = p;
 93	      }	
 94	      break;
 95	    }
 96	  }
 97	}
 98	if(matched) break;
 99      }
100      if(!matched) vetoEvent;
101      Particle piP,piM;
102      matched = false;
103      if (rhoP.children().size()!=2) vetoEvent;
104      if (rhoM.children().size()!=2) vetoEvent;
105      if(rhoP.children()[0].pid()==211 &&
106	 rhoP.children()[1].pid()==111)
107	piP = rhoP.children()[0];
108      else if(rhoP.children()[1].pid()==211 &&
109	      rhoP.children()[0].pid()==111)
110	piP = rhoP.children()[1];
111      else
112	vetoEvent;
113      if(rhoM.children()[0].pid()==-211 &&
114	 rhoM.children()[1].pid()==111)
115	piM = rhoM.children()[0];
116      else if(rhoM.children()[1].pid()==-211 &&
117	      rhoM.children()[0].pid()==111)
118	piM = rhoM.children()[1];
119      else
120	vetoEvent;
121      // boost to the rho+ rest frame
122      LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(rhoP.momentum().betaVec());
123      Vector3 e1z = rhoP.momentum().p3().unit();
124      Vector3 e1y = e1z.cross(axis).unit();
125      Vector3 e1x = e1y.cross(e1z).unit();
126      Vector3 axis1 = boost1.transform(piP.momentum()).p3().unit();
127      double n1x(e1x.dot(axis1)),n1y(e1y.dot(axis1)),n1z(e1z.dot(axis1));
128      _h[0]->fill(n1z);
129      _h[2]->fill(atan2(n1y,n1x));
130      // boost to the rho- frame
131      LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(rhoM.momentum().betaVec());
132      Vector3 axis2 = boost2.transform(piM.momentum()).p3().unit();
133      double n2x(e1x.dot(axis2)),n2y(e1y.dot(axis2)),n2z(e1z.dot(axis2));
134      _h[1]->fill(n2z);
135      _h[3]->fill(atan2(n2y,n2x));
136      _h[4]->fill(axis.dot(rhoP.momentum().p3().unit()));
137    }
138
139
140    /// Normalise histograms etc., after the run
141    void finalize() {
142      for(unsigned int ix=0;ix<5;++ix)
143	normalize(_h[ix],1.,false);
144    }
145
146    /// @}
147
148
149    /// @name Histograms
150    /// @{
151    Histo1DPtr _h[5];
152    /// @}
153
154
155  };
156
157
158  RIVET_DECLARE_PLUGIN(BABAR_2008_I789011);
159
160}