1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6#include "Rivet/Projections/MissingMomentum.hh"
  7#include "Rivet/Projections/PromptFinalState.hh"
  8#include "Rivet/Projections/UnstableParticles.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// @brief  Lambda_c -> Lambda l+ nu_l asymmetry
 14  class ARGUS_1994_I371613 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(ARGUS_1994_I371613);
 19
 20
 21    /// @name Analysis methods
 22    /// @{
 23
 24    /// Book histograms and initialise projections before the run
 25    void init() {
 26      // Initialise and register projections
 27      declare(UnstableParticles(), "UFS" );
 28
 29      // Book histograms
 30      book(_h_Lambda, "/TMP/hLambda", 20,-1.,1.);
 31
 32    }
 33
 34    void findChildren(Particle parent, int sign, unsigned int & npart,
 35		      Particles & lambda, Particles & e, Particles & nu) {
 36      for(const Particle & child : parent.children()) {
 37	if(child.pid()==sign*PID::LAMBDA) {
 38	  lambda.push_back(child);
 39	  ++npart;
 40	}
 41	else if(child.pid()==-sign*PID::EMINUS || child.pid()==-sign*PID::MUON) {
 42	  e.push_back(child);
 43	  ++npart;
 44	}
 45	else if(child.pid()==sign*PID::NU_E || child.pid()==sign*PID::NU_MU) {
 46	  nu.push_back(child);
 47	  ++npart;
 48	}
 49	else if(!child.children().empty()) {
 50	  findChildren(child,sign,npart,lambda,e,nu);
 51	}
 52	else {
 53	  ++npart;
 54	}
 55      }
 56    }
 57
 58
 59    /// Perform the per-event analysis
 60    void analyze(const Event& event) {
 61      // loop over Lambda_c baryons
 62      for( const Particle& Lambdac : apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==4122)) {
 63	int sign = Lambdac.pid()/4122;
 64        Particles lambda,e,nu;
 65        unsigned int npart(0);
 66        findChildren(Lambdac,sign,npart,lambda,e,nu);
 67	if(npart!=3 || lambda.size()!=1 || e.size()!=1 || nu.size()!=1) continue;
 68	Particle baryon2;
 69	if(lambda[0].children()[0].pid()== sign*2212 &&
 70	   lambda[0].children()[1].pid()== -sign*211) {
 71	  baryon2 = lambda[0].children()[0];
 72	}
 73	else if(lambda[0].children()[1].pid()== sign*2212 &&
 74		lambda[0].children()[0].pid()== -sign*211) {
 75	  baryon2 = lambda[0].children()[1];
 76	}
 77	else
 78	  continue;
 79	// mass cut
 80	double mLL = (lambda[0].momentum()+e[0].momentum()).mass();
 81	if(mLL<1.85 || mLL>2.2) continue;
 82	// first boost to the Lambdac rest frame
 83	LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(Lambdac.momentum().betaVec());
 84	FourMomentum pbaryon1 = boost1.transform(lambda[0].momentum());
 85	FourMomentum pbaryon2 = boost1.transform(baryon2  .momentum());
 86	// to lambda rest frame
 87	LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pbaryon1.betaVec());
 88	Vector3 axis = pbaryon1.p3().unit();
 89	FourMomentum pp = boost2.transform(pbaryon2);
 90	// calculate angle
 91	double cTheta = pp.p3().unit().dot(axis);
 92	_h_Lambda->fill(cTheta);
 93      }
 94    }
 95
 96    pair<double,double> calcAlpha(Histo1DPtr hist) {
 97      if(hist->numEntries()==0.) return make_pair(0.,0.);
 98      double sum1(0.),sum2(0.);
 99      for (const auto& bin : hist->bins() ) {
100        double Oi = bin.sumW();
101        if(Oi==0.) continue;
102        double ai = 0.5*(bin.xMax()-bin.xMin());
103        double bi = 0.5*ai*(bin.xMax()+bin.xMin());
104        double Ei = bin.errW();
105        sum1 += sqr(bi/Ei);
106        sum2 += bi/sqr(Ei)*(Oi-ai);
107      }
108      return make_pair(sum2/sum1,sqrt(1./sum1));
109    }
110
111    /// Normalise histograms etc., after the run
112    void finalize() {
113      //  asymmetry
114      normalize(_h_Lambda);
115      Estimate1DPtr _h_alpha;
116      book(_h_alpha,1,1,1);
117      pair<double,double> alpha = calcAlpha(_h_Lambda);
118      _h_alpha->bin(1).set(alpha.first, alpha.second);
119    }
120
121    /// @}
122
123
124    /// @name Histograms
125    /// @{
126    Histo1DPtr _h_Lambda;
127    /// @}
128
129
130  };
131
132
133  RIVET_DECLARE_PLUGIN(ARGUS_1994_I371613);
134
135}