1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4#include "Rivet/Projections/DecayedParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief B- -> Lambdac+ pbar pi-
 10  class BABAR_2008_I791879 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2008_I791879);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // Initialise and register projections
 23      UnstableParticles ufs = UnstableParticles(Cuts::abspid==PID::BPLUS);
 24      declare(ufs, "UFS");
 25      DecayedParticles BP(ufs);
 26      BP.addStable( 4122);
 27      BP.addStable(-4122);
 28      declare(BP, "BP");
 29      // histos
 30      book(_h[0],"TMP/hA",refData(1,1,1));
 31      book(_h[1],2,1,1);
 32    }
 33
 34
 35    /// Perform the per-event analysis
 36    void analyze(const Event& event) {
 37      static const map<PdgId,unsigned int> & mode   = { {-4122,1}, { 2212,1}, { 211,1} };
 38      static const map<PdgId,unsigned int> & modeCC = { { 4122,1}, {-2212,1}, {-211,1} };
 39      DecayedParticles BP = apply<DecayedParticles>(event, "BP");
 40      for(unsigned int ix=0;ix<BP.decaying().size();++ix) {
 41	// select right decay mode
 42	int sign = 1;
 43	if      (BP.decaying()[ix].pid()>0 && BP.modeMatches(ix,3,mode  )) sign = 1;
 44	else if (BP.decaying()[ix].pid()<0 && BP.modeMatches(ix,3,modeCC)) sign =-1;
 45	else continue;
 46	const Particle & lam = BP.decayProducts()[ix].at(-4122*sign)[0];
 47	const Particle & pp  = BP.decayProducts()[ix].at( 2212*sign)[0];
 48	const Particle & pip = BP.decayProducts()[ix].at(  211*sign)[0];
 49	_h[0]->fill((lam.momentum()+pp.momentum()).mass());
 50	if(BP.decaying()[ix].children().size()==2 &&
 51	   (BP.decaying()[ix].children()[0].abspid()==4112 ||
 52	    BP.decaying()[ix].children()[1].abspid()==4112)) {
 53	  LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(BP.decaying()[ix].momentum().betaVec());
 54	  FourMomentum pSigma  = boost1.transform(lam.momentum()+pip.momentum());
 55	  Vector3 axis = -pSigma.p3().unit();
 56	  FourMomentum pLambda = boost1.transform(lam.momentum());
 57	  LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pSigma.betaVec());
 58	  _h[1]->fill(axis.dot(boost2.transform(pLambda).p3().unit()));
 59	}
 60      }
 61    }
 62
 63
 64    /// Normalise histograms etc., after the run
 65    void finalize() {
 66      // phase space volumes in the different bins (PR has python program which computed these)
 67      vector<double> phsp = {0.1351399774095767,0.24107906667337267,0.3088204600383566,0.36206863378662535,
 68			     0.40659958403332336,0.4450070368053994,0.47872719863744506,0.5086545750109557,
 69			     0.5353906093802568,0.5593621540798246,0.5808845790744056,0.6001982286826937,
 70			     0.6174908272498764,0.6329119405055217,0.6465826880700479,0.6586024850910634,
 71			     0.6690538534165392,0.67800593732165,0.685517125504117,0.6916370414215524,
 72			     0.6964080775759626,0.6998665941923083,0.7020438666137178,0.7029668415116569,
 73			     0.7026587454114315,0.701139577428797,0.6984265098419529,0.6945342141075089,
 74			     0.6894751254590167,0.6832596558172149,0.6758963620528058,0.6673920744258324,
 75			     0.657751988084915,0.6469797186756903,0.6350773212151493,0.6220452692510539,
 76			     0.6078823887065968,0.5925857373802724,0.5761504163391165,0.5585692926468113};
 77      // normalize the histograms
 78      for(unsigned int ix=0;ix<2;++ix)
 79	normalize(_h[ix],1.,false);
 80      // convert first plot to scatter and normalize to phase space volume in bin
 81      // convert to scatter
 82      Estimate1DPtr tmp;
 83      book(tmp,1,1,1);
 84      barchart(_h[0],tmp);
 85      double step = 0.0344;
 86      // divide by phase space volume
 87      for(unsigned int ip=0;ip<tmp->numBins();++ip)
 88	tmp->bin(ip+1).scale(1./phsp[ip]/step);
 89    }
 90    /// @}
 91
 92
 93    /// @name Histograms
 94    /// @{
 95    Histo1DPtr _h[2];
 96    /// @}
 97
 98
 99  };
100
101
102  RIVET_DECLARE_PLUGIN(BABAR_2008_I791879);
103
104}