#ifndef LESTER_MTSTAR_EXPORT_H
#define LESTER_MTSTAR_EXPORT_H

/**
  
Copyright Christopher Lester.  Please contact Alan Barr or Christopher Lester if you wish to incorporate it in your code.  We will grant permission (!) but just want to have an idea of whether the code is un use.

The code evaluates a kinematic variable, mTStar, providing a maximal lower bound on the mass of a resonance R undergoing the decay decay:

R -> IJ

followed by           I -> SP

and                   J -> TQ

where at least one of I or J has the mass "mOnShellIntermediate", and where
the intermediate mass constraint
is enforced for at least one side of the event, where P and Q are invisible and only constrained by
the sum of their transverse componente (ptmiss).


If you use this code, please cite the paper where we introduced this variable in the context  
of H -> WW -> lnu lnu decays:

http://arxiv.org/abs/1108.3468

A negative return value suggests an internal error in the implementation and should be reported to the authors together with a simple test program reproducing the problem.


To evaluate mTStar call the function 

double Lester::mTBound(...)

defined below.  Note that to use this code you must edit FlatRandom.h to make it contain a routine that returns high-quality random numbers in the range [0,1).

*/

namespace Lester {

  // fwd dec
  double mTStar(double se, double sx, double sy, double sz,
		double te, double tx, double ty, double tz,
		double pmissx, double pmissy,
		double mOnShellIntermediate);


}

#include "FlatRandom.h"

#include <cmath>
#include <cassert>
#include <iostream>



namespace Lester {


 const double aTypicalScaleAgainstWhichStarAlgWasDeveloped = 26.;


const bool DEBUG2STAR = false;



/// This is an implementation of what is written up in Lab Book 8B.  Roughly page 20. -- stop press -- no it isn't -- this is an implementation of a version that relaxes the mass constraint on side 2.
/// Higgs -> W (1) WStar (2)
/// W (1) -> s(vis) + p(invis) // here we do apply a W mass constraint
/// WStar (2) -> t(vis) + q(invis) // no W mass contraint applied
double HelperFunctionHiggsWWStarMassLesterAtFixedKxKy_relaxedOnSide2Only(
	const double kx, const double ky,
	const double se, const double sx, const double sy, const double sz,
	const double te, const double tx, const double ty, const double tz,
        const double pmissx, const double pmissy,
	const double mw, 
	bool & wasSillyRef) {

wasSillyRef = true;

const double mwsq = mw*mw;
const double mssq = se*se-sx*sx-sy*sy-sz*sz;
const double mtsq = te*te-tx*tx-ty*ty-tz*tz;

const double pkx1 = pmissx + kx;
const double pky1 = pmissy + ky;
// const double pkx2 = pmissx - kx;
// const double pky2 = pmissy - ky;

const double star1 = mwsq - mssq + sx*pkx1+sy*pky1;
// const double star2 = mwsq - mtsq + tx*pkx2+ty*pky2;

const double pk1sq = pkx1*pkx1 + pky1*pky1;
// const double pk2sq = pkx2*pkx2 + pky2*pky2;


const double etssq = se*se-sz*sz;
const double ettsq = te*te-tz*tz;

bool good = true;

const double disc1 = star1*star1 - etssq*pk1sq;
// const double disc2 = star2*star2 - ettsq*pk2sq;

if (disc1<0) {
	good = false;
}


// where is invis 1?
const double px=0.5*(pmissx + kx);
const double py=0.5*(pmissy + ky);
const double pz_cat = (star1*sz + se*sqrt(disc1))/(2*etssq);
const double pz_dog = (star1*sz - se*sqrt(disc1))/(2*etssq);
const double pesq_cat = px*px + py*py + pz_cat*pz_cat;
const double pesq_dog = px*px + py*py + pz_dog*pz_dog;
const double pe_cat = sqrt(pesq_cat);
const double pe_dog = sqrt(pesq_dog);
 


// find 4-mom of everything else other than q:

const double alle_cat = pe_cat + se + te;
const double alle_dog = pe_dog + se + te;
const double allx     = px     + sx + tx;
const double ally     = py     + sy + ty;
const double allz_cat = pz_cat + sz + tz;
const double allz_dog = pz_dog + sz + tz;

// find mass of everything else other than q:

const double allmsq_cat = alle_cat*alle_cat - allx*allx - ally*ally - allz_cat*allz_cat;
const double allmsq_dog = alle_dog*alle_dog - allx*allx - ally*ally - allz_dog*allz_dog;

// find et of everything else other than q:

const double alletsq_cat = fabs(alle_cat*alle_cat - allz_cat*allz_cat); // have put the protection in, since alletsq_cat SHOULD be positive if the vectors have any pt.
const double alletsq_dog = fabs(alle_dog*alle_dog - allz_dog*allz_dog); // have put the protection in, since alletsq_dog SHOULD be positive if the vectors have any pt.
const double allet_cat = sqrt(alletsq_cat);
const double allet_dog = sqrt(alletsq_dog);

// find q_T
const double qx=0.5*(pmissx - kx);
const double qy=0.5*(pmissy - ky);
const double etq = sqrt(qx*qx + qy*qy);


// we know mass of q_T is zero by construction, so can now calculate the transverse mass of the all + q system treating "all" as visible (even though it contains a bit of p).
const double transversemasssq_cat = allmsq_cat + 0 + 2.0*( allet_cat*etq - allx*qx - ally*qy );
const double transversemasssq_dog = allmsq_dog + 0 + 2.0*( allet_dog*etq - allx*qx - ally*qy );

const double mhAAAsq = transversemasssq_cat;
const double mhBBBsq = transversemasssq_dog;

if (!good) {
        double power=1; // power to raise the -ve discriminant by in attemt to get out of bad region
	return pow(-disc1,power) + 14000;
}

const double mhsq = ((mhAAAsq <= mhBBBsq) ? mhAAAsq : mhBBBsq );

if (mhsq<0) {
	return -sqrt(-mhsq);
}

wasSillyRef = false;

return sqrt(mhsq);


} // end of function



/// This is an implementation of what is written up in Lab Book 8B.  Roughly page 20. -- stop press -- no it isn't -- this is an implementation of a version that relaxes the mass constraint on side 2.
// BOTH
/// Higgs -> W (1) WStar (2)
/// W (1) -> s(vis) + p(invis) // here we do apply a W mass constraint
/// WStar (2) -> t(vis) + q(invis) // no W mass contraint applied
// AND
/// Higgs -> W (1) WStar (2)
/// W (1) -> s(vis) + p(invis) // no W mass constraint applied
/// WStar (2) -> t(vis) + q(invis) // here we do apply a W mass contraint
// are considered, and the smaller is returned.


double HelperFunctionHiggsWWStarMassLesterAtFixedKxKy(
	const double kx, const double ky,
	const double se, const double sx, const double sy, const double sz,
	const double te, const double tx, const double ty, const double tz,
        const double pmissx, const double pmissy,
	const double mw, 
	bool & wasSillyRef) {

  wasSillyRef = true;

  bool wasSilly1;
  bool wasSilly2;

  const double mh1 = HelperFunctionHiggsWWStarMassLesterAtFixedKxKy_relaxedOnSide2Only(+kx,+ky,
									 se,sx,sy,sz,
									 te,tx,ty,tz,
									 pmissx,pmissy,
									 mw,
									 wasSilly1);
  const double mh2 = HelperFunctionHiggsWWStarMassLesterAtFixedKxKy_relaxedOnSide2Only(-kx,-ky,
									 te,tx,ty,tz,
									 se,sx,sy,sz,
									 pmissx,pmissy,
									 mw,
									 wasSilly2);

  if (wasSilly1 && wasSilly2) {
    return (mh1<=mh2) ? mh1 : mh2;
  }
  if (wasSilly1 && !wasSilly2) {
    wasSillyRef = false;
    return mh2;
  }
  if (wasSilly2 && !wasSilly1) {
    wasSillyRef = false;
    return mh1;
  }
  if (!wasSilly2 && !wasSilly1) {
    wasSillyRef = false;
    return (mh1<=mh2) ? mh1 : mh2;
  }
  return -20; // error code
}



  ////////////////////////////////////////////////////

// A negative return value indicates an error!
  double mTStar(double se, double sx, double sy, double sz,
		double te, double tx, double ty, double tz,
		double pmissx, double pmissy,
		double mOnShellIntermediate) {



    const double scale = sqrt((	sx*sx + sy*sy + sz*sz + 
				tx*tx + ty*ty + tz*tz + 
				pmissx*pmissx + pmissy*pmissy)/8.);
    // std::cout << "SCALE " << scale << std::endl;
    if (scale==0) {
      return mOnShellIntermediate;
    }


    // OK .. scale was not zero, so let's condition the inputs by dividing out this scale so that we are insensitive to unit choices that users make:

    se /= scale;
    sx /= scale;
    sy /= scale;
    sz /= scale;

    te /= scale;
    tx /= scale;
    ty /= scale;
    tz /= scale;

    pmissx /= scale;
    pmissy /= scale;

    mOnShellIntermediate /= scale;

    // We will have to multiply that scale back into our answer at the end!


    const double mw = mOnShellIntermediate;

    double kxStart = 0;
    double kyStart = 0;
    double bestHMassSoFar;

    bool haveValidStartPoint=false;
    bool bestPointWasSilly;

    // Attempt to get valid start point

    const double distFromWall /* AKA scan size */ = 2./aTypicalScaleAgainstWhichStarAlgWasDeveloped; // order of magnitude spread over which cauchy vals will be distributed.
    
    for (int i=0; i<10000; ++i) {

      // Warning! This doggedly tries to get a start point .. and may take forever if there is no possible start point ... or if the start point is 
	

      const double theta = (FlatRandom()-0.5)*3.14159;
      const double distToStep = distFromWall*tan(theta); 
      const double angToStep = FlatRandom()*3.14159*2.0;
      const double kx = distToStep * cos(angToStep);
      const double ky = distToStep * sin(angToStep);
      bool wasSilly;

      const double possHMass =
	Lester::HelperFunctionHiggsWWStarMassLesterAtFixedKxKy(kx,ky,
						 se,sx,sy,sz,
						 te,tx,ty,tz,
						 pmissx,pmissy,
						 mw,
						 wasSilly);
      if (   (possHMass>=0 && !haveValidStartPoint) ||  (possHMass>=0 && haveValidStartPoint && possHMass < bestHMassSoFar)   )  {
	bestHMassSoFar = possHMass;
	haveValidStartPoint = true;
	bestPointWasSilly = wasSilly;
	kxStart = kx;
	kyStart = ky;
	if (Lester::DEBUG2STAR) { 
	    std::cout << "Found better start point with mass " << possHMass << std::endl;
        }
	if (!wasSilly) {
	  // Don't need to work any harder ...
	  break;
	}
      }
	
     
    }
    if (Lester::DEBUG2STAR) {
	std::cout << "LESTER ANS " << bestHMassSoFar << "\tsilliness " << bestPointWasSilly << " valid " << haveValidStartPoint << std::endl;
    }


    if (haveValidStartPoint) {
      // Now we can attempt to minimise this function.
      
      double kxOld = kxStart;
      double kyOld = kyStart;
      double oldH = bestHMassSoFar;

      double typicalStepSize = distFromWall;
      const double shrinkageFactor = 0.99;
      const double growthFactor = 2;
      

      while(typicalStepSize > 1.e-6/aTypicalScaleAgainstWhichStarAlgWasDeveloped) {
	const double theta = (FlatRandom()-0.5)*3.14159;
	const double distToStep = typicalStepSize*tan(theta); 
	const double angToStep = FlatRandom()*3.14159*2.0001;
	const double newkx = kxOld + distToStep * cos(angToStep);
	const double newky = kyOld + distToStep * sin(angToStep);
	bool wasSilly;
 
	const double possHMass =
	  Lester::HelperFunctionHiggsWWStarMassLesterAtFixedKxKy(newkx,newky,
						   se,sx,sy,sz,
						   te,tx,ty,tz,
						   pmissx,pmissy,
						   mw,
						   wasSilly);
	if ( possHMass>=0 && possHMass < bestHMassSoFar) {
	  bestHMassSoFar = possHMass;
	  bestPointWasSilly = wasSilly;
	  kxOld = newkx;
	  kyOld = newky;
	  typicalStepSize *= growthFactor;
	  if (Lester::DEBUG2STAR) {
		std::cout << "Found better point with mass " << possHMass << " while typical step " << typicalStepSize << std::endl;
          }
	  
	} // if point is an improvement
	else {
	  typicalStepSize *= shrinkageFactor;
	} // point was not an improvement
	
      } // while we wantto keep going
      if (bestPointWasSilly) {
	return -10;
      } else {
	return bestHMassSoFar*scale;
      }
    } // have valid start
    else {
      return -15;
    } // have invalid start

 
 

}
 



} // end of namespace Lester

#endif