file ColliderBit/topness.h

[No description available]

Classes

Name
structmy_func
classmy_simplex
classmy_Nelder_Mead

Functions

Name
doublemy_lp(double u[4], double v[4])
voidmy_add(double a[4], double b[4], double c[4])
doublemy_dot(double p1[], double p2[], int d)
doublemy_dot(double p1[], double p2[], double a[], int d)
doublemy_enorm(double p[], int d)
doubletopnesscompute(double pb1[4], double pl[4], double MET[4], double sigmat, double sigmaW)

Attributes

Name
const doublemt
const doublemW
const intDMAX

Functions Documentation

function my_lp

double my_lp(
    double u[4],
    double v[4]
)

function my_add

void my_add(
    double a[4],
    double b[4],
    double c[4]
)

function my_dot

double my_dot(
    double p1[],
    double p2[],
    int d
)

function my_dot

double my_dot(
    double p1[],
    double p2[],
    double a[],
    int d
)

function my_enorm

double my_enorm(
    double p[],
    int d
)

function topnesscompute

double topnesscompute(
    double pb1[4],
    double pl[4],
    double MET[4],
    double sigmat,
    double sigmaW
)

Attributes Documentation

variable mt

const double mt =172.;

variable mW

const double mW =80.4;

variable DMAX

const int DMAX =4;

Source code

#include <iostream>
#include <string>
#include <cstdlib>
#include <cmath>

#include "gambit/Utils/threadsafe_rng.hpp"

// ****** Copy from https://github.com/michaelgraesser/topness *****
// ****** arXiv:1212.4495 ******************************************
// Modified according to https://arxiv.org/pdf/1612.03877.pdf
/* "The definition of topness used in this analysis is modified from
    the one originally proposed in [arXiv:1212.4495]: namely, the terms
    corresponding to the detected leptonic top quark decay and the
    centre-of-mass energy are dropped since in events with low jet
    multiplicity the second b jet is often not identified.
    In these cases, the discriminating power of the topness variable
    is reduced when a light-flavour jet is used instead in the calculation.

    Modified further by Pat Scott, Feb 9 2019, to use GAMBIT random
    number generator.

    This is a BSD License. Code written by Michael L. Graesser.

    Copyright (c) 2016, Los Alamos National Security, LLC
    All rights reserved.
    Copyright 2016. Los Alamos National Security, LLC. This software was produced under U.S. Government contract DE-AC52-06NA25396 for Los Alamos National Laboratory (LANL), which is operated by Los Alamos National Security, LLC for the U.S. Department of Energy. The U.S. Government has rights to use, reproduce, and distribute this software.  NEITHER THE GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.  If software is modified to produce derivative works, such modified software should be clearly marked, so as not to confuse it with the version available from LANL.

    Additionally, redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
    1.       Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
    2.      Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
    3.      Neither the name of Los Alamos National Security, LLC, Los Alamos National Laboratory, LANL, the U.S. Government, nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

    THIS SOFTWARE IS PROVIDED BY LOS ALAMOS NATIONAL SECURITY, LLC AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL LOS ALAMOS NATIONAL SECURITY, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */

// using namespace std;
// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
// >>>>>> topness_struct.cpp & topness_struct.h in "topness" <<<<<<<
const double mt=172.;  // top quark mass
const double mW=80.4;  // W mass
double my_lp(double u[4], double v[4])
{    // Lorentz product . assume x, y Lorentz vectors ordered as (px,py,pz,E) etc. Initialize with
  // metric g==(-1,-1,-1,1)

  double g[]={-1.,-1.,-1.,1.};
  double s=0;
  for (int j=0;j<4;j++)
    {s+= u[j]*g[j]*v[j]; }
  return s;
}
void my_add(double a[4], double b[4], double c[4])
{
  for (int i=0; i<4; i++)
    {
      c[i]=a[i]+b[i];
    }
}
struct my_func
{
  double pb1[4];
//  double pb2[4];
  double pl[4];
  double pMET[4];
  double sigmat, sigmaW, sign;
  // ***Modified by Yang Zhang 23.1.2018
  my_func(double ppb1[4], double ppl[4], double ppMET[4], double sigmatt, double sigmaWW, double ssign) :
  sigmat (sigmatt), sigmaW(sigmaWW), sign(ssign) {
    pb1[0]=ppb1[0];
    pb1[1]=ppb1[1];
    pb1[2]=ppb1[2];
    pb1[3]=ppb1[3];
//    pb2[0]=ppb2[0];
//    pb2[1]=ppb2[1];
//    pb2[2]=ppb2[2];
//    pb2[3]=ppb2[3];
    pl[0]=ppl[0];
    pl[1]=ppl[1];
    pl[2]=ppl[2];
    pl[3]=ppl[3];
    pMET[0]=ppMET[0];
    pMET[1]=ppMET[1];
    pMET[2]=ppMET[2];
    pMET[3]=ppMET[3];
  }
  double operator()(double points[],int /*d*/) {  // points[0]=pv_x, points[1]=pv_y, points[2]=pv_z, points[3]=pW_z
    // d is size of points = 4
    // pv_x, pv_y, pv_z
    double pvx=points[0];
    double pvy=points[1];
    double pvz=points[2];
    // neutrino energy assuming mass-shell condition
    double Ev=sqrt(pow(pvx,2)+pow(pvy,2)+pow(pvz,2));

    double pv[]={pvx,pvy,pvz,Ev};

    // pW_z
    double pWz=points[3];
    // std::cout << "points=" << pvx << ", " << pvy << ", " << pvz << ", " << pWz << std::endl;
    // W momenta from neutrino and MET
    double pW[]={-pvx+pMET[0],-pvy+pMET[1],pWz,sqrt(pow(-pvx+pMET[0],2)+pow(-pvy+pMET[1],2)+pow(pWz,2) + pow(mW,2))};

    double pb1W[4];
    my_add(pb1,pW,pb1W);
    double plv[4];
    my_add(pl,pv,plv);
//    double pb2lv[4];
//    my_add(pb2,plv,pb2lv);
    // function to minimize
    double fsum=0.;
    // ***Modified by Yang Zhang 23.1.2018
    fsum= pow(my_lp(pb1W,pb1W)-pow(mt,2),2)/pow(sigmat,4)+pow(my_lp(plv,plv)-pow(mW,2),2)/pow(sigmaW,4);
    //+pow(my_lp(pb2lv,pb2lv)-pow(mt,2),2)/pow(sigmat,4)+pow(my_lp(plv,plv)-pow(mW,2),2)/pow(sigmaW,4);
    //std::cout << "fsum = " << fsum << std::endl;
    return sign*fsum;
  }
};

// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
// >>>>>> simplex.cpp & simplex.h in "my_Nelder_Mead" <<<<<<<<<<<<<<
const int DMAX=4; // # parameters to scan over, not the space-time dimension! Number of points in the simplex is DMAX+1.
double my_dot(double p1[], double p2[], int d)   // d is size of arrays
{
// compute euclidean dot product between p1 and p2
  double value=0.;
  for (int i=0; i<d; i++)
    {
      //      std::cout << "p1: " << p1[i] <<", p2:" << p2[i] << std::endl;
      value+=p1[i]*p2[i];
    }
  //std::cout << "dp value is " << value << std::endl;
  return value;
}
double my_dot(double p1[], double p2[], double a[], int d)   // d is size of arrays
{
  // compute euclidean dot product between p1 and p2 with weights a

  double value=0.;
  for (int i=0; i<d; i++)
    {
      value+=a[i]*p1[i]*p2[i];
    }
  return value;
}
double my_enorm(double p[], int d)
{
  // compute euclidean norm of p
  double value=0.;
  double norm=0.;
  for (int i=0; i<d; i++)
    {
      value+=p[i]*p[i];
    }
  if (value>0.)
    {
      norm=sqrt(value);
    }
  else
    {
      norm=0.;
    }
  return norm;
}
class my_simplex{
// co-ordinate space is d dimensional. analysis uses N=d+1 points
 private:
  int d;
  double alpha, beta, gamma;

 public:
  my_func *f;
  my_simplex(int, double, double, double, my_func *);
  double xstart[DMAX*(DMAX+1)]; // length (d+1)*d. d elements are coordinates
  double x[DMAX*(DMAX+1)]; //    current polygon instance containing d+1 points of d dimension each, so size is d(d+1)
  double xh[DMAX]; // coordinates of point with highest value
  double xl[DMAX]; // coordinates of point with lowest value
  double y[DMAX+1]; // function values at x
  double yl, ynh, yh; // lowest, next highest, and highest function values
  double xCentroid[DMAX]; // current centroid;
  double yReflect, yExpand, yContract; // function values at these points
  double xReflect[DMAX]; // P* reflection
  double xExpand[DMAX]; // P** expansion
  double xContract[DMAX]; // P*** contraction
  void find_min();
  void find_max();
  int imin,imax, inmax;
  void my_SetUp(double xin[]);
  void set_y();
  void my_Centroid(int);
  void my_Reflection();
  void my_Expansion();
  void my_Contraction();
  void replace_all();
  double get_yavg();
  double get_sigma();
  void print_Centroid();
  void print_Reflect();
  void print_Expand();
  void print_Contract();
  void print_max();
  void print_min();
  void print_all();
  void print_xy();
  void print_xyl();
  void print_xyh();
};

class my_Nelder_Mead {
 private:
  int d;
  int Ntry;      // N cycles of Nelder-Mead algorithm
  double eps;

 public:
  my_func *f;
  my_simplex simplex;
  bool convergeYes; // set to True when algorithm has converged
  my_Nelder_Mead(int, double, double, double, int, double, my_func *);
  bool one_cycle(my_simplex *);
  bool find_global_min(double xin[DMAX*(DMAX+1)]);
  double yfinal;
  double xfinal[DMAX*(DMAX+1)];
};


my_simplex::my_simplex(int dd, double aalpha, double bbeta, double ggamma,  my_func (*ff)) : d(dd), alpha(aalpha), beta(bbeta), gamma(ggamma), f(ff){}
// check d and size of xin agree

void my_simplex::my_SetUp(double xin[DMAX*(DMAX+1)])
{
  int D=d*(d+1);
  std::copy(xin,xin+D,xstart);
  for (int i=0; i<D; i++)
    {
      x[i]=xstart[i];
    }

  std::vector<double> xi(d);
  for (int i=0; i<d+1; i++)
    {
      std::copy(x+d*i,x+d*i+d,&xi[0]);  // get coordinates of i'th point and copy into xi
      y[i]=(*f)(&xi[0],d);
    }
}

void my_simplex::set_y()
{
  std::vector<double> xi(d);
  for (int i=0; i<d+1; i++)
    {
      std::copy(x+d*i,x+d*i+d,&xi[0]);  // get coordinates of i'th point and copy into xi
      y[i]=(*f)(&xi[0],d);
    }
}

void my_simplex::find_max()
{
  for (int i=0; i<d+1; i++)
    {
      if (y[i]> y[imax])
    imax=i;
    }
  // now find second highest
  if (imax==1)
    inmax=0;

  for (int i=0; i<d+1; i++)
    {
      if (i==imax) continue;
      if (y[i] > y[inmax])
    {
      inmax=i;
    }
    }

  yh=y[imax];
  ynh=y[inmax];
  // find coordinates with maximum value
  for (int j=0; j<d; j++)
    {
      xh[j]=x[d*imax+j];
    }
}

void my_simplex::find_min()
{
  imin=0;   // first function value f(x_1)
  for (int i=0; i<d+1; i++)
    {
       if (y[i]< y[imin])
       {
        imin=i;
       }
    }

  yl=y[imin];

  // find coordinates with minimum value
  for (int j=0; j<d; j++)
    {
      xl[j]=x[d*imin+j];
    }

}

void my_simplex::my_Centroid(int h)
/* return P_bar */
{
  // compute centroid
  for (int j=0; j<d; j++)
    {
      xCentroid[j]=0;
      for (int i=0; i< d+1; i++)
    {
      if (i==h) continue;
      xCentroid[j]+=x[d*i+j];
    }
    }
}

void my_simplex::my_Reflection() // h is highest point
/* return P*
   REFLECTION  P*=(1+alpha)P_bar - alpha*P_h   */
{
  for (int i=0; i<d;i++)
  {
    xReflect[i] = xCentroid[i]*(1+alpha)/d- alpha*xh[i];
  }
  yReflect=(*f)(xReflect,d);
}

void my_simplex:: my_Expansion()
/* return EXPANSION P** =(1+gamma)*P* -gamma* P_bar */
{
  for (int j=0; j<d; j++)
  {
    xExpand[j] =xCentroid[j]*(1-gamma)/d+ (gamma)*xReflect[j];
  }
  yExpand=(*f)(xExpand,d);

}

void my_simplex:: my_Contraction()
/* return CONTRACTION P*** =beta* Ph +(1-beta)*P_bar */
{
  for (int j=0; j<d; j++)
  {
    xContract[j]=xCentroid[j]*(1-beta)/d+beta*xh[j];
  }
  yContract=(*f)(xContract,d);
}

void my_simplex::replace_all()
{
  for (int i=0; i<d+1; i++)
    {
      if (i==imin) continue;
      for (int j=0; j<d; j++)
    {
      x[d*i+j]=0.5*(x[d*i+j]+xl[j]);
    }
    }
}

double my_simplex::get_yavg()
{
  double yavg=0.;
  for (int i=0; i<d+1; i++)
    {
      yavg+=y[i];
    }
  yavg=yavg/(d+1);
  return yavg;
}

double my_simplex::get_sigma()
{
  double yavg=0.;
  double sigma=0;
  for (int i=0; i<d+1; i++)
    {
      yavg+=y[i];
    }
  yavg=yavg/(d+1);
  for (int i=0; i<d+1; i++)
    {
      sigma+=pow((y[i]-yavg),2);
    }

  sigma=sigma/(d+1);
  return sigma;
}

void my_simplex::print_Centroid()
{

  std::cout << "Current xCentroid is : " << std::endl;
  for (int k=0;k<d; k++)
    {
      std::cout << xCentroid[k];
      if (k==d-1)
    std::cout << std::endl;
      else std::cout<< ", ";

    }
}

void my_simplex::print_Reflect()
{
  std::cout << "Current xReflect and y value are : " << std::endl;
  for (int k=0;k<d; k++)
    {
      std::cout << xReflect[k];
      if (k==d-1)
    std::cout <<", " << yReflect << std::endl;
      else std::cout<< ", ";

    }
}

void my_simplex::print_Expand()
{

  std::cout << "Current xExpand and y value are : " << std::endl;
  for (int k=0;k<d; k++)
    {
      std::cout << xExpand[k];
      if (k==d-1)
    std::cout << ", " << yExpand << std::endl;
      else std::cout<< ", ";

    }
}

void my_simplex::print_Contract()
{
  std::cout << "Current xContract and y value are : " << std::endl;
  for (int k=0; k<d; k++)
    {
      std::cout << xContract[k];
      if (k==d-1)
    std::cout << ", " << yContract << std::endl;
      else std::cout << ", ";

    }

}

void my_simplex::print_max()
{
  std::cout << "Printing imax and inmax and their values " << std::endl;
  std::cout << "imax = " << imax << ", y[imax] = " << yh << std::endl;
  std::cout << "inmax = " << inmax << ", y[inmax] = " << ynh << std::endl;
}

void my_simplex::print_min()
{
  std::cout << "Printing imin and its value " << std::endl;
  std::cout << "imin = " << imin << ", y[imin] = " << yl << std::endl;
}

void my_simplex::print_xyh()
{
  std::cout << "The highest value is " << std::endl;
  for (int i=0; i<d; i++)
    {
      std::cout << xh[i];
      if ((i+1)% d !=0)
        {
          std::cout << ", " ;
        }
      else
        {
          std::cout << ", " << yh << std::endl;
        }
    }
}

void my_simplex::print_xyl()
{
  std::cout << "The lowest value is " << std::endl;
  for (int i=0; i<d; i++)
    {
      std::cout << xl[i];
      if ((i+1)% d !=0)
        {
          std::cout << ", " ;
        }
      else
        {
          std::cout << ", " << yl << std::endl;
        }
    }

}

void my_simplex::print_xy()
{
  // print current x and y
  std::cout << "Current x and y values are: " << std::endl;
  for (int i=0; i< d*(d+1); i++)
    {
      std::cout << x[i];
      if ((i+1) % d !=0)
        {
          std::cout << ", ";
        }
      else
        {
          div_t ratio;
          ratio=div(i,d);
      std::cout << ", " << y[ratio.quot] << std::endl;
        }
    }

}
void my_simplex::print_all()
{
// print current x
  std::cout << "Current x values are: " << std::endl;
  for (int i=0; i< d*(d+1); i++)
    {
      std::cout << x[i];
      if ((i+1) % d !=0)
    {
      std::cout << ", ";
    }
      else
    {
      div_t ratio;
      ratio=div(i,d);
      std::cout << ", " << y[ratio.quot]<<std::endl;
    }
    }
  std::cout << "Current centroid: " << std::endl;
  for (int i=0; i< d; i++)
      {
    std::cout << xCentroid[i];
          if ((i+1) % d !=0)
          {
          std::cout << ", ";
           }
          else
          {
        std::cout << std::endl;
          }
     }
  std::cout << "Current Reflection: " << std::endl;
  for (int i=0; i< d; i++)
     {
         std::cout << xReflect[i];
     if ((i+1) % d !=0)
     {
         std::cout << ", ";
     }
     else
     {
       std::cout << ", " << yReflect << std::endl;
     }
    }
  std::cout << "Current Expansion: " << std::endl;
  for (int i=0; i< d; i++)
      {
    std::cout << xExpand[i];
      if ((i+1) % d !=0)
      {
          std::cout << ", ";
      }
      else
      {
        std::cout << ", " << yExpand << std::endl;
      }
      }
  std::cout << "Current Contraction: " << std::endl;
  for (int i=0; i< d; i++)
      {
    std::cout << xContract[i];
      if ((i+1) % d !=0)
      {
          std::cout << ", ";
      }
      else
      {
        std::cout << ", " << yContract << std::endl;
      }
     }
  std::cout << "The highest value is " << std::endl;
  for (int i=0; i<d; i++)
    {
      std::cout << xh[i];
      if ((i+1)% d !=0)
    {
      std::cout << ", " ;
    }
      else
    {
      std::cout << ", " << yh << std::endl;
    }
    }
  std::cout << "The lowest value is " << std::endl;
  for (int i=0; i<d; i++)
    {
      std::cout << xl[i];
      if ((i+1)% d !=0)
        {
          std::cout << ", " ;
        }
      else
        {
          std::cout << ", " << yl << std::endl;
        }
    }
    std::cout << "The lowest point is imin=" << imin << std::endl;

    std::cout << "The highest point is imax=" << imax << std::endl;

    std::cout << "The next highest point is inmax=" << inmax << std::endl;
}


my_Nelder_Mead::my_Nelder_Mead(int dd, double alpha, double beta, double gamma, int NNtry, double eeps, my_func *ff): d(dd), Ntry(NNtry), eps(eeps), f(ff), simplex(dd, alpha, beta, gamma, f){}

bool my_Nelder_Mead::one_cycle(my_simplex *s)
{
// execute one-iteration of Nelder-Mead method
  (*s).my_Centroid((*s).imax);
  //  (*s).print_Centroid();
  (*s).my_Reflection();
  // (*s).print_Reflect();
  if ((*s).yReflect <= (*s).yl)
    {
      // do expansion
      (*s).my_Expansion();
      //      (*s).print_Expand();
      if ((*s).yExpand < (*s).yl)
    {
      // replace P_h with P_**
      std::copy((*s).xExpand,(*s).xExpand+d, (*s).x+d*((*s).imax));
      (*s).set_y();
      return false;
    }
      else
    {
      // replace P_h with P_*
      std::copy((*s).xReflect,(*s).xReflect+d, (*s).x+d*((*s).imax));
      (*s).set_y();
      return true;
    }
    }
  else if (((*s).yReflect) >= (*s).ynh ) // do contraction
    {
      if (((*s).yReflect) < (*s).yh)
    {
      std::copy((*s).xReflect,(*s).xReflect+d, (*s).x+d*((*s).imax));
      (*s).set_y();
      (*s).find_max();
    }
      (*s).my_Contraction();
      //    (*s).print_Contract();
      if ((*s).yContract < (*s).yh)
    {
      std::copy((*s).xContract,(*s).xContract+d, (*s).x+d*((*s).imax));
      (*s).set_y();
      return false;

    }
      else
    {
      (*s).replace_all();
      (*s).set_y();
      return false;
    }
     }
  else {
    // replace P_h with P_*
    std::copy((*s).xReflect,(*s).xReflect+d, (*s).x+d*((*s).imax));
    (*s).set_y();
    return true;
  }


}

bool my_Nelder_Mead::find_global_min(double xin[DMAX*(DMAX+1)])
/* try Nelder-Mead cycle Ntry times; if values converge then restart using point near new
minimum  */
{
  // initialize
  yfinal=10000000000000.; // a very large number
  simplex.imax=0;
  simplex.inmax=1;
  simplex.imin=2;
  simplex.my_SetUp(xin);
  simplex.find_max();
  simplex.find_min();
  //  simplex.print_xy();
  //simplex.print_xyh();
  //simplex.print_xyl();
  //simplex.print_max();
  //simplex.print_min();
  convergeYes=false;
  // bool reflectYes=false;
  //  double xnew[DMAX*(DMAX+1)];
  for (int i=0; i<Ntry; i++)
    {
      //reflectYes=
      one_cycle(&simplex);
      //if (reflectYes==true) --i;
      simplex.find_max();
      simplex.find_min();
      double ynewmin=simplex.yl;
      double ynewmax=simplex.yh;
      //     std::cout << "i=" << i << std::endl;
      //std::cout << std::endl << std::endl << std::endl;
      //simplex.print_xy();
      //simplex.print_xyh();
      //simplex.print_xyl();
      //simplex.print_max();
      //simplex.print_min();
      //std::cout << std::endl << std::endl << std::endl;
      // double yavg=simplex.get_yavg();
      //  double sigma=simplex.get_sigma();
      if (std::abs(ynewmax -ynewmin)/(std::abs(ynewmax)+std::abs(ynewmin)+eps) < eps)
      {
        convergeYes=true;

        // save old xi, generate new xstarti = xbest/rndm + rndm*step
        if (ynewmin< yfinal)
          {
        std::copy(simplex.x,simplex.x+d*(d+1),xfinal);
        yfinal=ynewmin;
        //      std::cout << "final i=" << i << ", ymin = " << ynewmin << ", ymax = " << ynewmax << std::endl;
        break;
          }
    }
    }

  return convergeYes;
}



// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
// >>>>>> Wrappertopness.cpp & Wrappertopness.h in "topness" <<<<<<<
double topnesscompute(double pb1[4], double pl[4], double MET[4], double sigmat, double sigmaW)
{

  double alpha=1.0; // parameter for reflection
  double beta=0.5;  // parameter for contraction
  double gamma=2.0;  // parameter for expansion

  const int d=4;    // number of parameters to scan over, not the space-time dimension!
  const int DIMMAX=d*(d+1);  // dimension of simplex
  double xin[d+1][d];    // starting point
  double xstart[DIMMAX]; // algorithm stores d+1 points of simplex in a single array


  double eps=0.000002;   // tolerance
  double Deltastep=20.;  // initial spacing of points, in GeV
  int Ntry=100000;    // maximum number of Nelder-Mead cycles to perform for a given initial seed
  int Nattempts=15; //   number of initial starts

  double edir[]={1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1};
  double ybest=1000000000.; // a big number
  double ybest1=100000000.;
  double xbest1[]={1000.,1000.,1000.,1000.};
//  double ybest2=100000000.;  // a big number
//  double xbest2[]={1000.,1000.,1000.,1000.};
  int i,j,k;
  // int tid;
  bool converge1;
// converge2;

  // initialize topness function  ***Modified by Yang Zhang 23.1.2018
  my_func topstat1(pb1,pl,MET,sigmat,sigmaW,1.0);
//  my_func topstat2(pb2,pb1,pl,MET,sigmat,sigmaW,1.0); // other combination

  // initialize topness computation
  my_Nelder_Mead my_check1(d,alpha,beta,gamma, Ntry, eps, &topstat1);
//  my_Nelder_Mead my_check2(d,alpha,beta,gamma, Ntry, eps, &topstat2);
  double yl;
  // begin loop over Nattempts
  for (k=0; k<Nattempts; k++)
    {
      // initialize starting point
      for (i=0;i<d+1; i++)
        {
          for (j=0; j<d; j++)
            {
              if (i==0)
                {
                  xin[0][j]=8000.0*(Gambit::Random::draw()-0.5);
                }
              else
                {
                  xin[i][j]=xin[0][j]+Deltastep*edir[d*(i-1)+j];
                }
            }
          std::copy(xin[i],xin[i]+d, xstart+d*i);  // copy initial data into xstart
        }
      // now first combination
      converge1=my_check1.find_global_min(xstart);
      if (converge1==true)
        {

          yl=my_check1.yfinal;
          if (yl < ybest1)
            {
              ybest1=yl;
              std::copy(my_check1.xfinal,my_check1.xfinal+d,xbest1);
            }
        }
      else
        {
          std::cout << " Minimum not found...exiting " << std::endl;
        }
//      // now do second combination
//      converge2=my_check2.find_global_min(xstart);
//      if (converge2==true)
//        {
//          yl=my_check2.yfinal;
//          if (yl < ybest2)
//            {
//              ybest2=yl;
//              std::copy(my_check2.xfinal,my_check2.xfinal+d,xbest2);
//            }
//        }
//      else
//      {
//          std::cout << " Minimum not found...exiting " << std::endl;
//        }

    }

//  if (ybest1 < ybest2)
//    {
    ybest=ybest1;
//    std::copy(xbest1,xbest1+d,xbest);
//    }
//   else
//   {
//     ybest=ybest2;
//     std::copy(xbest2,xbest2+d,xbest);
//   }

  return ybest;

}

Updated on 2024-07-18 at 13:53:34 +0000