gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
topness.h
Go to the documentation of this file.
1 #include <iostream>
2 #include <string>
3 #include <cstdlib>
4 #include <cmath>
5 
7 
8 // ****** Copy from https://github.com/michaelgraesser/topness *****
9 // ****** arXiv:1212.4495 ******************************************
10 // Modified according to https://arxiv.org/pdf/1612.03877.pdf
11 /* "The definition of topness used in this analysis is modified from
12  the one originally proposed in [arXiv:1212.4495]: namely, the terms
13  corresponding to the detected leptonic top quark decay and the
14  centre-of-mass energy are dropped since in events with low jet
15  multiplicity the second b jet is often not identified.
16  In these cases, the discriminating power of the topness variable
17  is reduced when a light-flavour jet is used instead in the calculation.
18 
19  Modified further by Pat Scott, Feb 9 2019, to use GAMBIT random
20  number generator.
21 
22  This is a BSD License. Code written by Michael L. Graesser.
23 
24  Copyright (c) 2016, Los Alamos National Security, LLC
25  All rights reserved.
26  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.
27 
28  Additionally, redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
29  1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
30  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.
31  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.
32 
33  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.
34  */
35 
36 // using namespace std;
37 // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
38 // >>>>>> topness_struct.cpp & topness_struct.h in "topness" <<<<<<<
39 const double mt=172.; // top quark mass
40 const double mW=80.4; // W mass
41 double my_lp(double u[4], double v[4])
42 { // Lorentz product . assume x, y Lorentz vectors ordered as (px,py,pz,E) etc. Initialize with
43  // metric g==(-1,-1,-1,1)
44 
45  double g[]={-1.,-1.,-1.,1.};
46  double s=0;
47  for (int j=0;j<4;j++)
48  {s+= u[j]*g[j]*v[j]; }
49  return s;
50 }
51 void my_add(double a[4], double b[4], double c[4])
52 {
53  for (int i=0; i<4; i++)
54  {
55  c[i]=a[i]+b[i];
56  }
57 }
58 struct my_func
59 {
60  double pb1[4];
61 // double pb2[4];
62  double pl[4];
63  double pMET[4];
64  double sigmat, sigmaW, sign;
65  // ***Modified by Yang Zhang 23.1.2018
66  my_func(double ppb1[4], double ppl[4], double ppMET[4], double sigmatt, double sigmaWW, double ssign) :
67  sigmat (sigmatt), sigmaW(sigmaWW), sign(ssign) {
68  pb1[0]=ppb1[0];
69  pb1[1]=ppb1[1];
70  pb1[2]=ppb1[2];
71  pb1[3]=ppb1[3];
72 // pb2[0]=ppb2[0];
73 // pb2[1]=ppb2[1];
74 // pb2[2]=ppb2[2];
75 // pb2[3]=ppb2[3];
76  pl[0]=ppl[0];
77  pl[1]=ppl[1];
78  pl[2]=ppl[2];
79  pl[3]=ppl[3];
80  pMET[0]=ppMET[0];
81  pMET[1]=ppMET[1];
82  pMET[2]=ppMET[2];
83  pMET[3]=ppMET[3];
84  }
85  double operator()(double points[],int /*d*/) { // points[0]=pv_x, points[1]=pv_y, points[2]=pv_z, points[3]=pW_z
86  // d is size of points = 4
87  // pv_x, pv_y, pv_z
88  double pvx=points[0];
89  double pvy=points[1];
90  double pvz=points[2];
91  // neutrino energy assuming mass-shell condition
92  double Ev=sqrt(pow(pvx,2)+pow(pvy,2)+pow(pvz,2));
93 
94  double pv[]={pvx,pvy,pvz,Ev};
95 
96  // pW_z
97  double pWz=points[3];
98  // std::cout << "points=" << pvx << ", " << pvy << ", " << pvz << ", " << pWz << std::endl;
99  // W momenta from neutrino and MET
100  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))};
101 
102  double pb1W[4];
103  my_add(pb1,pW,pb1W);
104  double plv[4];
105  my_add(pl,pv,plv);
106 // double pb2lv[4];
107 // my_add(pb2,plv,pb2lv);
108  // function to minimize
109  double fsum=0.;
110  // ***Modified by Yang Zhang 23.1.2018
111  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);
112  //+pow(my_lp(pb2lv,pb2lv)-pow(mt,2),2)/pow(sigmat,4)+pow(my_lp(plv,plv)-pow(mW,2),2)/pow(sigmaW,4);
113  //std::cout << "fsum = " << fsum << std::endl;
114  return sign*fsum;
115  }
116 };
117 
118 // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
119 // >>>>>> simplex.cpp & simplex.h in "my_Nelder_Mead" <<<<<<<<<<<<<<
120 const int DMAX=4; // # parameters to scan over, not the space-time dimension! Number of points in the simplex is DMAX+1.
121 double my_dot(double p1[], double p2[], int d) // d is size of arrays
122 {
123 // compute euclidean dot product between p1 and p2
124  double value=0.;
125  for (int i=0; i<d; i++)
126  {
127  // std::cout << "p1: " << p1[i] <<", p2:" << p2[i] << std::endl;
128  value+=p1[i]*p2[i];
129  }
130  //std::cout << "dp value is " << value << std::endl;
131  return value;
132 }
133 double my_dot(double p1[], double p2[], double a[], int d) // d is size of arrays
134 {
135  // compute euclidean dot product between p1 and p2 with weights a
136 
137  double value=0.;
138  for (int i=0; i<d; i++)
139  {
140  value+=a[i]*p1[i]*p2[i];
141  }
142  return value;
143 }
144 double my_enorm(double p[], int d)
145 {
146  // compute euclidean norm of p
147  double value=0.;
148  double norm=0.;
149  for (int i=0; i<d; i++)
150  {
151  value+=p[i]*p[i];
152  }
153  if (value>0.)
154  {
155  norm=sqrt(value);
156  }
157  else
158  {
159  norm=0.;
160  }
161  return norm;
162 }
164 // co-ordinate space is d dimensional. analysis uses N=d+1 points
165  private:
166  int d;
167  double alpha, beta, gamma;
168 
169  public:
171  my_simplex(int, double, double, double, my_func *);
172  double xstart[DMAX*(DMAX+1)]; // length (d+1)*d. d elements are coordinates
173  double x[DMAX*(DMAX+1)]; // current polygon instance containing d+1 points of d dimension each, so size is d(d+1)
174  double xh[DMAX]; // coordinates of point with highest value
175  double xl[DMAX]; // coordinates of point with lowest value
176  double y[DMAX+1]; // function values at x
177  double yl, ynh, yh; // lowest, next highest, and highest function values
178  double xCentroid[DMAX]; // current centroid;
179  double yReflect, yExpand, yContract; // function values at these points
180  double xReflect[DMAX]; // P* reflection
181  double xExpand[DMAX]; // P** expansion
182  double xContract[DMAX]; // P*** contraction
183  void find_min();
184  void find_max();
186  void my_SetUp(double xin[]);
187  void set_y();
188  void my_Centroid(int);
189  void my_Reflection();
190  void my_Expansion();
191  void my_Contraction();
192  void replace_all();
193  double get_yavg();
194  double get_sigma();
195  void print_Centroid();
196  void print_Reflect();
197  void print_Expand();
198  void print_Contract();
199  void print_max();
200  void print_min();
201  void print_all();
202  void print_xy();
203  void print_xyl();
204  void print_xyh();
205 };
206 
208  private:
209  int d;
210  int Ntry; // N cycles of Nelder-Mead algorithm
211  double eps;
212 
213  public:
216  bool convergeYes; // set to True when algorithm has converged
217  my_Nelder_Mead(int, double, double, double, int, double, my_func *);
218  bool one_cycle(my_simplex *);
219  bool find_global_min(double xin[DMAX*(DMAX+1)]);
220  double yfinal;
221  double xfinal[DMAX*(DMAX+1)];
222 };
223 
224 
225 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){}
226 // check d and size of xin agree
227 
228 void my_simplex::my_SetUp(double xin[DMAX*(DMAX+1)])
229 {
230  int D=d*(d+1);
231  std::copy(xin,xin+D,xstart);
232  for (int i=0; i<D; i++)
233  {
234  x[i]=xstart[i];
235  }
236 
237  for (int i=0; i<d+1; i++)
238  {
239  double xi[d];
240  std::copy(x+d*i,x+d*i+d,xi); // get coordinates of i'th point and copy into xi
241  y[i]=(*f)(xi,d);
242  }
243 }
244 
246 {
247  for (int i=0; i<d+1; i++)
248  {
249  double xi[d];
250  std::copy(x+d*i,x+d*i+d,xi); // get coordinates of i'th point and copy into xi
251  y[i]=(*f)(xi,d);
252  }
253 }
254 
256 {
257  for (int i=0; i<d+1; i++)
258  {
259  if (y[i]> y[imax])
260  imax=i;
261  }
262  // now find second highest
263  if (imax==1)
264  inmax=0;
265 
266  for (int i=0; i<d+1; i++)
267  {
268  if (i==imax) continue;
269  if (y[i] > y[inmax])
270  {
271  inmax=i;
272  }
273  }
274 
275  yh=y[imax];
276  ynh=y[inmax];
277  // find coordinates with maximum value
278  for (int j=0; j<d; j++)
279  {
280  xh[j]=x[d*imax+j];
281  }
282 }
283 
285 {
286  imin=0; // first function value f(x_1)
287  for (int i=0; i<d+1; i++)
288  {
289  if (y[i]< y[imin])
290  {
291  imin=i;
292  }
293  }
294 
295  yl=y[imin];
296 
297  // find coordinates with minimum value
298  for (int j=0; j<d; j++)
299  {
300  xl[j]=x[d*imin+j];
301  }
302 
303 }
304 
306 /* return P_bar */
307 {
308  // compute centroid
309  for (int j=0; j<d; j++)
310  {
311  xCentroid[j]=0;
312  for (int i=0; i< d+1; i++)
313  {
314  if (i==h) continue;
315  xCentroid[j]+=x[d*i+j];
316  }
317  }
318 }
319 
320 void my_simplex::my_Reflection() // h is highest point
321 /* return P*
322  REFLECTION P*=(1+alpha)P_bar - alpha*P_h */
323 {
324  for (int i=0; i<d;i++)
325  {
326  xReflect[i] = xCentroid[i]*(1+alpha)/d- alpha*xh[i];
327  }
328  yReflect=(*f)(xReflect,d);
329 }
330 
332 /* return EXPANSION P** =(1+gamma)*P* -gamma* P_bar */
333 {
334  for (int j=0; j<d; j++)
335  {
336  xExpand[j] =xCentroid[j]*(1-gamma)/d+ (gamma)*xReflect[j];
337  }
338  yExpand=(*f)(xExpand,d);
339 
340 }
341 
343 /* return CONTRACTION P*** =beta* Ph +(1-beta)*P_bar */
344 {
345  for (int j=0; j<d; j++)
346  {
347  xContract[j]=xCentroid[j]*(1-beta)/d+beta*xh[j];
348  }
349  yContract=(*f)(xContract,d);
350 }
351 
353 {
354  for (int i=0; i<d+1; i++)
355  {
356  if (i==imin) continue;
357  for (int j=0; j<d; j++)
358  {
359  x[d*i+j]=0.5*(x[d*i+j]+xl[j]);
360  }
361  }
362 }
363 
365 {
366  double yavg=0.;
367  for (int i=0; i<d+1; i++)
368  {
369  yavg+=y[i];
370  }
371  yavg=yavg/(d+1);
372  return yavg;
373 }
374 
376 {
377  double yavg=0.;
378  double sigma=0;
379  for (int i=0; i<d+1; i++)
380  {
381  yavg+=y[i];
382  }
383  yavg=yavg/(d+1);
384  for (int i=0; i<d+1; i++)
385  {
386  sigma+=pow((y[i]-yavg),2);
387  }
388 
389  sigma=sigma/(d+1);
390  return sigma;
391 }
392 
394 {
395 
396  std::cout << "Current xCentroid is : " << std::endl;
397  for (int k=0;k<d; k++)
398  {
399  std::cout << xCentroid[k];
400  if (k==d-1)
401  std::cout << std::endl;
402  else std::cout<< ", ";
403 
404  }
405 }
406 
408 {
409  std::cout << "Current xReflect and y value are : " << std::endl;
410  for (int k=0;k<d; k++)
411  {
412  std::cout << xReflect[k];
413  if (k==d-1)
414  std::cout <<", " << yReflect << std::endl;
415  else std::cout<< ", ";
416 
417  }
418 }
419 
421 {
422 
423  std::cout << "Current xExpand and y value are : " << std::endl;
424  for (int k=0;k<d; k++)
425  {
426  std::cout << xExpand[k];
427  if (k==d-1)
428  std::cout << ", " << yExpand << std::endl;
429  else std::cout<< ", ";
430 
431  }
432 }
433 
435 {
436  std::cout << "Current xContract and y value are : " << std::endl;
437  for (int k=0; k<d; k++)
438  {
439  std::cout << xContract[k];
440  if (k==d-1)
441  std::cout << ", " << yContract << std::endl;
442  else std::cout << ", ";
443 
444  }
445 
446 }
447 
449 {
450  std::cout << "Printing imax and inmax and their values " << std::endl;
451  std::cout << "imax = " << imax << ", y[imax] = " << yh << std::endl;
452  std::cout << "inmax = " << inmax << ", y[inmax] = " << ynh << std::endl;
453 }
454 
456 {
457  std::cout << "Printing imin and its value " << std::endl;
458  std::cout << "imin = " << imin << ", y[imin] = " << yl << std::endl;
459 }
460 
462 {
463  std::cout << "The highest value is " << std::endl;
464  for (int i=0; i<d; i++)
465  {
466  std::cout << xh[i];
467  if ((i+1)% d !=0)
468  {
469  std::cout << ", " ;
470  }
471  else
472  {
473  std::cout << ", " << yh << std::endl;
474  }
475  }
476 }
477 
479 {
480  std::cout << "The lowest value is " << std::endl;
481  for (int i=0; i<d; i++)
482  {
483  std::cout << xl[i];
484  if ((i+1)% d !=0)
485  {
486  std::cout << ", " ;
487  }
488  else
489  {
490  std::cout << ", " << yl << std::endl;
491  }
492  }
493 
494 }
495 
497 {
498  // print current x and y
499  std::cout << "Current x and y values are: " << std::endl;
500  for (int i=0; i< d*(d+1); i++)
501  {
502  std::cout << x[i];
503  if ((i+1) % d !=0)
504  {
505  std::cout << ", ";
506  }
507  else
508  {
509  div_t ratio;
510  ratio=div(i,d);
511  std::cout << ", " << y[ratio.quot] << std::endl;
512  }
513  }
514 
515 }
517 {
518 // print current x
519  std::cout << "Current x values are: " << std::endl;
520  for (int i=0; i< d*(d+1); i++)
521  {
522  std::cout << x[i];
523  if ((i+1) % d !=0)
524  {
525  std::cout << ", ";
526  }
527  else
528  {
529  div_t ratio;
530  ratio=div(i,d);
531  std::cout << ", " << y[ratio.quot]<<std::endl;
532  }
533  }
534  std::cout << "Current centroid: " << std::endl;
535  for (int i=0; i< d; i++)
536  {
537  std::cout << xCentroid[i];
538  if ((i+1) % d !=0)
539  {
540  std::cout << ", ";
541  }
542  else
543  {
544  std::cout << std::endl;
545  }
546  }
547  std::cout << "Current Reflection: " << std::endl;
548  for (int i=0; i< d; i++)
549  {
550  std::cout << xReflect[i];
551  if ((i+1) % d !=0)
552  {
553  std::cout << ", ";
554  }
555  else
556  {
557  std::cout << ", " << yReflect << std::endl;
558  }
559  }
560  std::cout << "Current Expansion: " << std::endl;
561  for (int i=0; i< d; i++)
562  {
563  std::cout << xExpand[i];
564  if ((i+1) % d !=0)
565  {
566  std::cout << ", ";
567  }
568  else
569  {
570  std::cout << ", " << yExpand << std::endl;
571  }
572  }
573  std::cout << "Current Contraction: " << std::endl;
574  for (int i=0; i< d; i++)
575  {
576  std::cout << xContract[i];
577  if ((i+1) % d !=0)
578  {
579  std::cout << ", ";
580  }
581  else
582  {
583  std::cout << ", " << yContract << std::endl;
584  }
585  }
586  std::cout << "The highest value is " << std::endl;
587  for (int i=0; i<d; i++)
588  {
589  std::cout << xh[i];
590  if ((i+1)% d !=0)
591  {
592  std::cout << ", " ;
593  }
594  else
595  {
596  std::cout << ", " << yh << std::endl;
597  }
598  }
599  std::cout << "The lowest value is " << std::endl;
600  for (int i=0; i<d; i++)
601  {
602  std::cout << xl[i];
603  if ((i+1)% d !=0)
604  {
605  std::cout << ", " ;
606  }
607  else
608  {
609  std::cout << ", " << yl << std::endl;
610  }
611  }
612  std::cout << "The lowest point is imin=" << imin << std::endl;
613 
614  std::cout << "The highest point is imax=" << imax << std::endl;
615 
616  std::cout << "The next highest point is inmax=" << inmax << std::endl;
617 }
618 
619 
620 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){}
621 
623 {
624 // execute one-iteration of Nelder-Mead method
625  (*s).my_Centroid((*s).imax);
626  // (*s).print_Centroid();
627  (*s).my_Reflection();
628  // (*s).print_Reflect();
629  if ((*s).yReflect <= (*s).yl)
630  {
631  // do expansion
632  (*s).my_Expansion();
633  // (*s).print_Expand();
634  if ((*s).yExpand < (*s).yl)
635  {
636  // replace P_h with P_**
637  std::copy((*s).xExpand,(*s).xExpand+d, (*s).x+d*((*s).imax));
638  (*s).set_y();
639  return false;
640  }
641  else
642  {
643  // replace P_h with P_*
644  std::copy((*s).xReflect,(*s).xReflect+d, (*s).x+d*((*s).imax));
645  (*s).set_y();
646  return true;
647  }
648  }
649  else if (((*s).yReflect) >= (*s).ynh ) // do contraction
650  {
651  if (((*s).yReflect) < (*s).yh)
652  {
653  std::copy((*s).xReflect,(*s).xReflect+d, (*s).x+d*((*s).imax));
654  (*s).set_y();
655  (*s).find_max();
656  }
657  (*s).my_Contraction();
658  // (*s).print_Contract();
659  if ((*s).yContract < (*s).yh)
660  {
661  std::copy((*s).xContract,(*s).xContract+d, (*s).x+d*((*s).imax));
662  (*s).set_y();
663  return false;
664 
665  }
666  else
667  {
668  (*s).replace_all();
669  (*s).set_y();
670  return false;
671  }
672  }
673  else {
674  // replace P_h with P_*
675  std::copy((*s).xReflect,(*s).xReflect+d, (*s).x+d*((*s).imax));
676  (*s).set_y();
677  return true;
678  }
679 
680 
681 }
682 
684 /* try Nelder-Mead cycle Ntry times; if values converge then restart using point near new
685 minimum */
686 {
687  // initialize
688  yfinal=10000000000000.; // a very large number
689  simplex.imax=0;
690  simplex.inmax=1;
691  simplex.imin=2;
692  simplex.my_SetUp(xin);
693  simplex.find_max();
694  simplex.find_min();
695  // simplex.print_xy();
696  //simplex.print_xyh();
697  //simplex.print_xyl();
698  //simplex.print_max();
699  //simplex.print_min();
700  convergeYes=false;
701  // bool reflectYes=false;
702  // double xnew[DMAX*(DMAX+1)];
703  for (int i=0; i<Ntry; i++)
704  {
705  //reflectYes=
706  one_cycle(&simplex);
707  //if (reflectYes==true) --i;
708  simplex.find_max();
709  simplex.find_min();
710  double ynewmin=simplex.yl;
711  double ynewmax=simplex.yh;
712  // std::cout << "i=" << i << std::endl;
713  //std::cout << std::endl << std::endl << std::endl;
714  //simplex.print_xy();
715  //simplex.print_xyh();
716  //simplex.print_xyl();
717  //simplex.print_max();
718  //simplex.print_min();
719  //std::cout << std::endl << std::endl << std::endl;
720  // double yavg=simplex.get_yavg();
721  // double sigma=simplex.get_sigma();
722  if (std::abs(ynewmax -ynewmin)/(std::abs(ynewmax)+std::abs(ynewmin)+eps) < eps)
723  {
724  convergeYes=true;
725 
726  // save old xi, generate new xstarti = xbest/rndm + rndm*step
727  if (ynewmin< yfinal)
728  {
729  std::copy(simplex.x,simplex.x+d*(d+1),xfinal);
730  yfinal=ynewmin;
731  // std::cout << "final i=" << i << ", ymin = " << ynewmin << ", ymax = " << ynewmax << std::endl;
732  break;
733  }
734  }
735  }
736 
737  return convergeYes;
738 }
739 
740 
741 
742 // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
743 // >>>>>> Wrappertopness.cpp & Wrappertopness.h in "topness" <<<<<<<
744 double topnesscompute(double pb1[4], double pl[4], double MET[4], double sigmat, double sigmaW)
745 {
746 
747  double alpha=1.0; // parameter for reflection
748  double beta=0.5; // parameter for contraction
749  double gamma=2.0; // parameter for expansion
750 
751  const int d=4; // number of parameters to scan over, not the space-time dimension!
752  const int DIMMAX=d*(d+1); // dimension of simplex
753  double xin[d+1][d]; // starting point
754  double xstart[DIMMAX]; // algorithm stores d+1 points of simplex in a single array
755 
756 
757  double eps=0.000002; // tolerance
758  double Deltastep=20.; // initial spacing of points, in GeV
759  int Ntry=100000; // maximum number of Nelder-Mead cycles to perform for a given initial seed
760  int Nattempts=15; // number of initial starts
761 
762  double edir[]={1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1};
763  double ybest=1000000000.; // a big number
764  double ybest1=100000000.;
765  double xbest1[]={1000.,1000.,1000.,1000.};
766 // double ybest2=100000000.; // a big number
767 // double xbest2[]={1000.,1000.,1000.,1000.};
768  int i,j,k;
769  // int tid;
770  bool converge1;
771 // converge2;
772 
773  // initialize topness function ***Modified by Yang Zhang 23.1.2018
774  my_func topstat1(pb1,pl,MET,sigmat,sigmaW,1.0);
775 // my_func topstat2(pb2,pb1,pl,MET,sigmat,sigmaW,1.0); // other combination
776 
777  // initialize topness computation
778  my_Nelder_Mead my_check1(d,alpha,beta,gamma, Ntry, eps, &topstat1);
779 // my_Nelder_Mead my_check2(d,alpha,beta,gamma, Ntry, eps, &topstat2);
780  double yl;
781  // begin loop over Nattempts
782  for (k=0; k<Nattempts; k++)
783  {
784  // initialize starting point
785  for (i=0;i<d+1; i++)
786  {
787  for (j=0; j<d; j++)
788  {
789  if (i==0)
790  {
791  xin[0][j]=8000.0*(Gambit::Random::draw()-0.5);
792  }
793  else
794  {
795  xin[i][j]=xin[0][j]+Deltastep*edir[d*(i-1)+j];
796  }
797  }
798  std::copy(xin[i],xin[i]+d, xstart+d*i); // copy initial data into xstart
799  }
800  // now first combination
801  converge1=my_check1.find_global_min(xstart);
802  if (converge1==true)
803  {
804 
805  yl=my_check1.yfinal;
806  if (yl < ybest1)
807  {
808  ybest1=yl;
809  std::copy(my_check1.xfinal,my_check1.xfinal+d,xbest1);
810  }
811  }
812  else
813  {
814  std::cout << " Minimum not found...exiting " << std::endl;
815  }
816 // // now do second combination
817 // converge2=my_check2.find_global_min(xstart);
818 // if (converge2==true)
819 // {
820 // yl=my_check2.yfinal;
821 // if (yl < ybest2)
822 // {
823 // ybest2=yl;
824 // std::copy(my_check2.xfinal,my_check2.xfinal+d,xbest2);
825 // }
826 // }
827 // else
828 // {
829 // std::cout << " Minimum not found...exiting " << std::endl;
830 // }
831 
832  }
833 
834 // if (ybest1 < ybest2)
835 // {
836  ybest=ybest1;
837 // std::copy(xbest1,xbest1+d,xbest);
838 // }
839 // else
840 // {
841 // ybest=ybest2;
842 // std::copy(xbest2,xbest2+d,xbest);
843 // }
844 
845  return ybest;
846 
847 }
bool convergeYes
Definition: topness.h:216
void print_max()
Definition: topness.h:448
void replace_all()
Definition: topness.h:352
double x[DMAX *(DMAX+1)]
Definition: topness.h:173
double pl[4]
Definition: topness.h:62
START_MODEL beta
Definition: Axions.hpp:36
int imin
Definition: topness.h:185
bool find_global_min(double xin[DMAX *(DMAX+1)])
Definition: topness.h:683
double my_enorm(double p[], int d)
Definition: topness.h:144
my_func(double ppb1[4], double ppl[4], double ppMET[4], double sigmatt, double sigmaWW, double ssign)
Definition: topness.h:66
double yfinal
Definition: topness.h:220
double yExpand
Definition: topness.h:179
void print_min()
Definition: topness.h:455
double xCentroid[DMAX]
Definition: topness.h:178
double pb1[4]
Definition: topness.h:60
const double mt
Definition: topness.h:39
void find_min()
Definition: topness.h:284
double get_yavg()
Definition: topness.h:364
double sign
Definition: topness.h:64
double yl
Definition: topness.h:177
void print_all()
Definition: topness.h:516
double pMET[4]
Definition: topness.h:63
void find_max()
Definition: topness.h:255
double get_sigma()
Definition: topness.h:375
double xExpand[DMAX]
Definition: topness.h:181
START_MODEL b
Definition: demo.hpp:270
const double mW
Definition: topness.h:40
my_func * f
Definition: topness.h:170
void print_xy()
Definition: topness.h:496
double xfinal[DMAX *(DMAX+1)]
Definition: topness.h:221
double yContract
Definition: topness.h:179
double sigmaW
Definition: topness.h:64
void print_Centroid()
Definition: topness.h:393
double yh
Definition: topness.h:177
double operator()(double points[], int)
Definition: topness.h:85
my_func * f
Definition: topness.h:214
double yReflect
Definition: topness.h:179
START_MODEL alpha
void my_Expansion()
Definition: topness.h:331
double y[DMAX+1]
Definition: topness.h:176
void print_Contract()
Definition: topness.h:434
double xl[DMAX]
Definition: topness.h:175
double sigmat
Definition: topness.h:64
static double draw()
Draw a single uniform random deviate from the interval (0,1) using the chosen RNG engine...
const double sigma
Definition: SM_Z.hpp:43
double xh[DMAX]
Definition: topness.h:174
void print_xyl()
Definition: topness.h:478
int imax
Definition: topness.h:185
my_simplex(int, double, double, double, my_func *)
Definition: topness.h:225
void my_Reflection()
Definition: topness.h:320
void print_Expand()
Definition: topness.h:420
double my_dot(double p1[], double p2[], int d)
Definition: topness.h:121
double ynh
Definition: topness.h:177
double alpha
Definition: topness.h:167
double eps
Definition: topness.h:211
int inmax
Definition: topness.h:185
A threadsafe interface to the STL random number generators.
double my_lp(double u[4], double v[4])
Definition: topness.h:41
void print_xyh()
Definition: topness.h:461
my_simplex simplex
Definition: topness.h:215
void print_Reflect()
Definition: topness.h:407
void my_SetUp(double xin[])
Definition: topness.h:228
double xReflect[DMAX]
Definition: topness.h:180
double pow(const double &a)
Outputs a^i.
void set_y()
Definition: topness.h:245
double xContract[DMAX]
Definition: topness.h:182
bool one_cycle(my_simplex *)
Definition: topness.h:622
double topnesscompute(double pb1[4], double pl[4], double MET[4], double sigmat, double sigmaW)
Definition: topness.h:744
double xstart[DMAX *(DMAX+1)]
Definition: topness.h:172
const int DMAX
Definition: topness.h:120
void my_add(double a[4], double b[4], double c[4])
Definition: topness.h:51
double beta
Definition: topness.h:167
void my_Centroid(int)
Definition: topness.h:305
double gamma
Definition: topness.h:167
void my_Contraction()
Definition: topness.h:342
my_Nelder_Mead(int, double, double, double, int, double, my_func *)
Definition: topness.h:620