1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief Excited Lambda_c spectra
 10  class ARGUS_1997_I440304 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ARGUS_1997_I440304);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // projections
 23      declare(Beam(), "Beams");
 24      declare(UnstableParticles(), "UFS");
 25      // book histos
 26      book(_h_rate1,"TMP/tot_1",refData(1,1,1));
 27      book(_h_rate2,"TMP/tot_2",refData(1,2,1));
 28      book(_h_x,2,1,1);
 29    }
 30
 31
 32    void findDecayProducts(Particle parent, Particles & Lambda_c, Particles & pions,unsigned int & nstable) {
 33      for(const Particle & p : parent.children()) {
 34	if(p.abspid()==4122) {
 35	  Lambda_c.push_back(p);
 36	  ++nstable;
 37	}
 38	else if(p.abspid()==PID::PIPLUS) {
 39	  pions.push_back(p);
 40	  ++nstable;
 41	}
 42	else if(!p.children().empty())
 43	  findDecayProducts(p,Lambda_c,pions,nstable);
 44	else
 45	  ++nstable;
 46      }
 47    }
 48
 49    /// Perform the per-event analysis
 50    void analyze(const Event& event) {
 51      static const int id2595 = 102142;
 52      // Get beams and average beam momentum
 53      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 54      const double Emax = ( beams.first.p3().mod() + beams.second.p3().mod() ) / 2.0;
 55      const double Pmax = sqrt(sqr(Emax)-sqr(2.595));
 56      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 57      for (const Particle& p : ufs.particles(Cuts::abspid==id2595)) {
 58	// spectrum
 59	double xp = p.momentum().p3().mod()/Pmax;
 60	_h_x->fill(xp);
 61	Particles Lambda_c,pions;
 62	unsigned int nstable(0);
 63	findDecayProducts(p,Lambda_c,pions,nstable);
 64	if(nstable==3&&pions.size()==2&&Lambda_c.size()==1) {
 65	  _h_rate1->fill(xp);
 66	  _h_rate2->fill(xp);
 67	}
 68      }
 69    }
 70
 71
 72    /// Normalise histograms etc., after the run
 73    void finalize() {
 74      normalize(_h_x);
 75      // br for lambda_c mode from pdg 2018
 76      double br = 0.0623;
 77      scale(_h_rate1,br*crossSection()/sumOfWeights()/picobarn);
 78      scale(_h_rate2,br*crossSection()/sumOfWeights()/picobarn);
 79      Estimate1DPtr tmp;
 80      book(tmp,1,1,1);
 81      barchart(_h_rate1,tmp);
 82      book(tmp,1,2,1);
 83      barchart(_h_rate2,tmp);
 84    }
 85
 86    /// @}
 87
 88
 89    /// @name Histograms
 90    /// @{
 91    Histo1DPtr _h_x,_h_rate1,_h_rate2;
 92    /// @}
 93
 94
 95  };
 96
 97
 98  RIVET_DECLARE_PLUGIN(ARGUS_1997_I440304);
 99
100}