// -*- C++ -*-
//
// HiggsMTmScale.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the HiggsMTmScale class.
//

#include "HiggsMTmScale.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"

using namespace Herwig;

HiggsMTmScale::HiggsMTmScale(): ShowerScale(0) {}

HiggsMTmScale::~HiggsMTmScale() {}

IBPtr HiggsMTmScale::clone() const {
  return new_ptr(*this);
}

IBPtr HiggsMTmScale::fullclone() const {
  return new_ptr(*this);
}

Energy2 HiggsMTmScale::renormalizationScale() const {

  Energy mHiggs = ZERO;
  Energy ptHiggs = ZERO;
  bool higgsFound = false;

  auto iData = mePartonData().begin() + 2;
  auto iMom = meMomenta().begin() + 2;
  for ( ; iData != mePartonData().end(); ++iData, ++iMom ) {
    if ( (*iData)->id() == 25 ) {
      higgsFound = true;
      mHiggs = (*iData)->mass();
      ptHiggs = iMom->perp();
    }
  }

  if ( !higgsFound )
    throw Exception() << "HiggsMTmScale::renormalizationScale(): No Higgs could be found. Check your setup."
		      << Exception::runerror;

  return 0.5*mHiggs*sqrt(sqr(mHiggs)/4.+sqr(ptHiggs));
  
}

Energy2 HiggsMTmScale::factorizationScale() const {
  return renormalizationScale();
}





// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
Energy2 HiggsMTmScale::showerScale() const {

  /*
       FACTORIZATION SCALE
  */
  if (ShowerScale == 0)
    {
     	return renormalizationScale();
    }
    
  /* 
     AVERAGE PT 
  */
  else if (ShowerScale ==1) {

    Energy ht=ZERO;
    auto iData = mePartonData().begin() + 2;
    auto iMom = meMomenta().begin() + 2;
    int numColoured=0;
    for ( ; iData != mePartonData().end(); ++iData, ++iMom ) {
      // pT of the coloured partons
      if(abs((*iData)->id()) <= 6 || (*iData)->id() == 21){
	numColoured++;
	ht+=iMom->perp();
      }
    }
    // Add a fake minimum, this is normally excluded when imposing cuts
    return sqr( max(ht/numColoured, 5._GeV));
  }

  /* 
     MIN RECONSTRUCTED JET PT
  */
  else if (ShowerScale ==2) {
    tcPDVector pd (mePartonData().begin() + 2, mePartonData().end());
    vector<LorentzMomentum> p (meMomenta().begin() + 2, meMomenta().end());
    tcPDPtr t1 = mePartonData()[0];
    tcPDPtr t2 = mePartonData()[1];
    tcCutsPtr cuts = lastCutsPtr();
    theJetFinder->cluster(pd, p, cuts, t1, t2);
    Energy2 minpt2 = ZERO;
    tcPDVector::const_iterator itpd = pd.begin();
    for (vector<LorentzMomentum>::const_iterator itp = p.begin() ;
	 itp != p.end(); ++itp, ++itpd )
      if ( theJetFinder->unresolvedMatcher()->check(**itpd) ) {
	if(minpt2 == ZERO) {
	  minpt2 = (*itp).perp2();
	}else{
	  minpt2 = min(minpt2,(*itp).perp2());
	}
      }
    // Add a fake minimum, this is normally excluded when imposing cuts
    return max(minpt2, sqr(5._GeV));
  }

  /* Geometric average of the vector bosons virtualities */
  /* This value of lamba5 corresponds to the one for PDF4LHC15_nnlo_30*/
  else {
    LorentzMomentum p1(ZERO,ZERO,ZERO,ZERO);
    LorentzMomentum p2(ZERO,ZERO,ZERO,ZERO);
    tcPDVector pd (mePartonData().begin(), mePartonData().end());
    vector<LorentzMomentum> p (meMomenta().begin(), meMomenta().end());
    tcPDVector::const_iterator pdata = pd.begin();
    vector<LorentzMomentum>::const_iterator mom = p.begin();
    int count=0;

    for ( ; mom != p.end(); ++pdata, ++mom ) {
      if((*pdata)->id() != 25){
	count ++;
	double sign =-1.;
	if(count <=2) sign=+1;
	if(mom->z() > ZERO){
	  p1 += sign*(*mom);
	}else{
	  p2 += sign*(*mom);
	}
      }    
    }

    Energy2 q1sq, q2sq, lamba5MSBsq;

    q1sq = max(sqr(p1.x())+sqr(p1.y())+sqr(p1.z())-sqr(p1.t()), sqr(1_GeV));
    q2sq = max(sqr(p2.x())+sqr(p2.y())+sqr(p2.z())-sqr(p2.t()), sqr(1_GeV));
    lamba5MSBsq = sqr(0.20836823616052982_GeV);

    return lamba5MSBsq * exp( sqrt( log(q1sq/lamba5MSBsq)* log(q2sq/lamba5MSBsq)));
  }
  
  
}


void HiggsMTmScale::persistentOutput(PersistentOStream & os) const {
  os << theJetFinder << ShowerScale;
}

void HiggsMTmScale::persistentInput(PersistentIStream & is, int) {
  is >> theJetFinder >> ShowerScale;
  
}


// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<HiggsMTmScale,MatchboxScaleChoice>
  describeHerwigHiggsMTmScale("Herwig::HiggsMTmScale", "HiggsMTmScale.so");

void HiggsMTmScale::Init() {

  static ClassDocumentation<HiggsMTmScale> documentation
    ("HiggsMTmScale implements scale choices related to 1/2 ht_parton.");

  
  static Reference<HiggsMTmScale,JetFinder> interfaceJetFinder
    ("JetFinder",
     "A reference to the jet finder.",
     &HiggsMTmScale::theJetFinder, false, false, true, false, false);

  
  static Switch<HiggsMTmScale, unsigned int> ifaceShowerScale
    ("ShowerScale",
     "Choice of the ShowerScale",
     &HiggsMTmScale::ShowerScale, 0, false, false);
  static SwitchOption fact
    (ifaceShowerScale,"fact","factorization scale", 0);
  static SwitchOption avpt
    (ifaceShowerScale,"avpt","average partons pT", 1);
  static SwitchOption minpt
    (ifaceShowerScale,"minpt","min reconstructed jet pT", 2);
  static SwitchOption virt
    (ifaceShowerScale,"virt","geometric average of the vector bosons pT", 3);


}