SusyScan/PlotScript/PlotTools.cc

#include "PlotTools.h"
#include "SusyScan.h"
#include "GeneratorMasses.h"

#include <cmath>
#include <algorithm>
#include <iostream>

#include "TGraph.h"
#include "TF1.h"
#include "TH2.h"
#include "TH2F.h"
#include "TObjArray.h"
#include "TPad.h"
#include "TCanvas.h"
#include "TRint.h"
#include "TROOT.h"

template<class T>
TGraph * PlotTools<T>::Line( double(*x)(const T*), double(*y)(const T*), 
                             double(*func)(const T*), const double mass, const double diff )
{
  TGraph * result = new TGraph(1);
  std::sort(_scan->begin(),_scan->end(),sort_by(x));
  int i=0;
  for (typename std::vector<T*>::const_iterator it=_scan->begin();it!=_scan->end();++it){
     if ( fabs(func( *it)-mass)<diff )
       result->SetPoint(i++, x(*it), y(*it));
  } 
  result->SetLineWidth(0.5);
  result->SetLineColor(kGray);
  return result;
}


template<class T>
void PlotTools<T>::Area( TH2*h, double(*x)(const T*), double(*y)(const T*), 
                         double(*f)(const T*) )
{
  //std::cout<< _scan->size() <<",  &scan="<<_scan<<std::endl;
  for (typename std::vector<T*>::const_iterator it=_scan->begin();it!=_scan->end();++it){
    h->SetBinContent( h->GetXaxis()->FindBin(x(*it)), 
                      h->GetYaxis()->FindBin(y(*it)), f(*it) );
                      
//   std::cout <<"m0="<< (*it)->Mzero<<", m1/2="<< (*it)->Mhalf<< ", sq="<<(*it)->MUL << ", gl="        <<(*it)->MGL <<std::endl;             
                      
  } 
}

template<class T>
void PlotTools<T>::Graph( TGraph*g, double(*x)(const T*), double(*y)(const T*),double ymin)
{
  unsigned i = g->GetN();
  std::sort(_scan->begin(),_scan->end(),sort_by(x));
  for (typename std::vector<T*>::const_iterator it=_scan->begin();it!=_scan->end();++it){
    if (y(*it)>=ymin) g->SetPoint(i++, x(*it), y(*it));
    //std::cout << i << ": x="<<x(*it)<< ", y="<<y(*it)<< std::endl;
  } 
}

template<class T>
std::vector<TGraph*> PlotTools<T>::GetContours005(TH2*h, int ncont)
{
   double conts[1];
   conts[0] = 0.05;
   if (!h) return std::vector<TGraph*>();
   TH2 * plot = (TH2*)h->Clone();
   plot->SetContour(1, conts);
   plot->SetFillColor(1);
   plot->Draw("CONT Z List");
   gPad->Update();
   TObjArray *contours = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
   int ncontours      = contours->GetSize();
   std::vector<TGraph*> result;
   for (int i=0;i<ncontours;++i){
     TList *list = (TList*)contours->At(i);
     TGraph* curv = (TGraph*)list->First();
     if (curv) result.push_back( curv );
     for(int j = 0; j < list->GetSize(); j++){
         curv = (TGraph*)list->After(curv); // Get Next graph
         if (curv) result.push_back( curv );
     }
   }  
   delete plot;
   std::sort(result.begin(),result.end(),sort_TGraph());
   return result;
}

template<class T>
TGraph * PlotTools<T>::GetContour005(TH2*h, int ncont, int flag)
{
   return (TGraph*)GetContours005(h, ncont).at(flag)->Clone();
}

template<class T>
std::vector<TGraph*> PlotTools<T>::GetContours(TH2*h, int ncont)
{
   if (!h) return std::vector<TGraph*>();
   TH2 * plot = (TH2*)h->Clone();
   plot->SetContour(ncont);
   plot->SetFillColor(1);
   plot->Draw("CONT Z List");
   gPad->Update();
   TObjArray *contours = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
   int ncontours      = contours->GetSize();
   std::vector<TGraph*> result;
   for (int i=0;i<ncontours;++i){
     TList *list = (TList*)contours->At(i);
     TGraph* curv = (TGraph*)list->First();
     if (curv) result.push_back( curv );
     for(int j = 0; j < list->GetSize(); j++){
         curv = (TGraph*)list->After(curv); // Get Next graph
         if (curv) result.push_back( curv );
     }
   }  
   delete plot;
   std::sort(result.begin(),result.end(),sort_TGraph());
   return result;
}

template<class T>
TGraph * PlotTools<T>::GetContour(TH2*h, int ncont, int flag)
{
   return (TGraph*)GetContours(h, ncont).at(flag)->Clone();
}

template<class T>
TGraph * PlotTools<T>::GetContour(TH2*h,double(*x)(const T*), double(*y)(const T*), 
                      double(*func)(const T*), int ncont, int flag,
                      int color, int style){
  TH2*hist = (TH2*)h->Clone();
  Area(hist, x, y, func);
  TGraph * graph = GetContour(hist, ncont, flag);
  graph->SetLineColor(color);
  graph->SetLineStyle(style);
  return graph;
}
 
template<class T>
TGraph * PlotTools<T>::GetContour005(TH2*h,double(*x)(const T*), double(*y)(const T*), 
                      double(*func)(const T*), int ncont, int flag,
                      int color, int style){
  TH2*hist = (TH2*)h->Clone();
  Area(hist, x, y, func);
  TGraph * graph = GetContour005(hist, ncont, flag);
  graph->SetLineColor(color);
  graph->SetLineStyle(style);
  return graph;
}
 
template<class T>
void PlotTools<T>::Print(double(*f)(const T*), double(*x)(const T*), double(*y)(const T*), TGraph*g, double p)
{
   for (typename std::vector<T*>::const_iterator it=_scan->begin();it!=_scan->end();++it){
     for (int j=0; j<g->GetN(); ++j) {
       double gx, gy;
       g->GetPoint(j, gx, gy);
       if ( (x(*it)-gx)*(x(*it)-gx) +  (y(*it)-gy)*(y(*it)-gy) < p*p)
         std::cout << x(*it) << ", " << y(*it) 
                   << " :: " << f(*it) << std::endl;
     }
   } 

}
template<class T>
void PlotTools<T>::Print(double(*f)(const T*), 
                         double(*x1)(const T*), double(*x2)(const T*), 
                         double(*x3)(const T*), double(*x4)(const T*), 
                         TGraph*g, double p)
{
   for (typename std::vector<T*>::const_iterator it=_scan->begin();it!=_scan->end();++it){
     for (int j=0; j<g->GetN(); ++j) {
       double gx, gy;
       g->GetPoint(j, gx, gy);
       if ( (x1(*it)-gx)*(x1(*it)-gx) +  (x2(*it)-gy)*(x2(*it)-gy) < p*p)
         std::cout << x1(*it) << ", " << x2(*it) << ", "  
                   << x3(*it) << ", " << x4(*it)
                   << " :: " << f(*it) << std::endl;
     }
   } 

}

template<class T>
TGraph * PlotTools<T>::ChooseBest(TGraph*r1,TGraph*r2,TGraph*g1,TGraph*g2,double x,double y)
{
  TGraph * res = new TGraph(0);
  for (int i=0; i<r1->GetN(); ++i) {
    double rx, r1y;
    r1->GetPoint(i,rx,r1y);
    double r2y = r2->Eval(rx);
    double g1y=g1->Eval(rx);
    double resy = (r1y>=r2y?g1y:g2->Eval(rx));
    if (rx>x&&r1y>y)
      res->SetPoint(i,rx,resy);
    else
      res->SetPoint(i,rx,g1y);
  }
  res->SetLineColor(g1->GetLineColor());
  res->SetLineWidth(g1->GetLineWidth());
  res->SetLineStyle(g1->GetLineStyle());
  res->SetFillColor(g1->GetFillColor());
  res->SetFillStyle(g1->GetFillStyle());
  return res;
}

template<class T>
TH2 * PlotTools<T>::BinWiseOr(TH2*h1, TH2*h2)
{
  TH2 * res = (TH2*)h1->Clone();
  for (int x=0; x<=res->GetXaxis()->GetNbins(); ++x)
    for (int y=0; y<=res->GetYaxis()->GetNbins(); ++y)
      if (h2->GetBinContent(x,y)>0.5)
        res->SetBinContent(x,y,h2->GetBinContent(x,y));
  return res;   
}

template<class T>
TGraph * PlotTools<T>::ModifyExpSigma(TGraph*g1, TGraph*g2, TGraph*g3)
{
  TGraph * res = new TGraph(0);
  for (int i=0; i<g1->GetN(); ++i) {
    double x, y;
    g1->GetPoint(i,x,y);
    double y2 = g2->Eval(x);
    double y3 = g3->Eval(x);
    res->SetPoint(i,x,y-y2+y3);
  }
  res->SetLineColor(g1->GetLineColor());
  res->SetLineStyle(g1->GetLineStyle());
  res->SetLineWidth(g1->GetLineWidth());
  res->SetFillColor(g1->GetFillColor());
  res->SetFillStyle(g1->GetFillStyle());
  return res;
}

double GetX(TGraph* g, double x0, double y0){
  double dn=999999999, up=9999999999, minx, maxx;
  for (int i=0; i<g->GetN(); ++i) {
    double x, y;
    g->GetPoint(i,x,y);
    if (y==y0) return x;
    if ( y<y0 && (x-x0)*(x-x0)+(y-y0)*(y-y0)<dn ){
      dn=(x-x0)*(x-x0)+(y-y0)*(y-y0);
      minx=x;
    }
    if ( y>y0 && (x-x0)*(x-x0)+(y-y0)*(y-y0)<up ){
      up=(x-x0)*(x-x0)+(y-y0)*(y-y0);
      maxx=x;
    }
  }
  return (minx+maxx)/2.;
}

template<class T>
TGraph * PlotTools<T>::ModifyExpSigmaY(TGraph*g1, TGraph*g2, TGraph*g3)
{
  TGraph * res = new TGraph(0);
  for (int i=0; i<g1->GetN(); ++i) {
    double x, y;
    g1->GetPoint(i,x,y);
    double x2,y2,x3,y3;
    x2=GetX(g2,x,y);
    x3=GetX(g3,x,y);
    res->SetPoint(i,x-x2+x3,y);
  }
  res->SetLineColor(g1->GetLineColor());
  res->SetLineStyle(g1->GetLineStyle());
  res->SetLineWidth(g1->GetLineWidth());
  res->SetFillColor(g1->GetFillColor());
  res->SetFillStyle(g1->GetFillStyle());
  return res;
}


template<class T>
bool PlotTools<T>::sort_TGraph::operator()(const TGraph*g1, const TGraph*g2)
{ 
   return g1->GetN() > g2->GetN();
}

TGraph * MakeBand(TGraph *g1, TGraph *g2, bool b){
  TGraph * res = new TGraph(g1->GetN()+g2->GetN());
  int p=0;
  for (int i=0; i<g1->GetN(); ++i) {
    double x, y;
    g1->GetPoint(i, x, y);
    res->SetPoint(p++, x, y);
  }
  for (int i=g2->GetN()-1; i>=0; --i) {
    double x, y;
    g2->GetPoint(i, x, y);
    res->SetPoint(p++, x, y);
  }
  res->SetLineColor( g1->GetLineColor() );
  res->SetFillColor( g2->GetLineColor() );
  res->SetFillStyle(4050);
  return res;
}

void Smooth(TGraph * g, int N)
{
  TGraph * old = (TGraph*)g->Clone();
  //int N = (n%2==0?n+1:n);
  if (N>2*g->GetN()) N=2*g->GetN()-1;


  double gauss[N];
  double sigma = (double)N/4.;
  double sum=0;
  double lim=(double)N/2.;
  TF1 *fb = new TF1("fb","gaus(0)",-lim,lim);
  fb->SetParameter(0, 1./(sqrt(2*3.1415)*sigma) );
  fb->SetParameter(1, 0);
  fb->SetParameter(2, sigma);
  for (int i=0; i<N; ++i){
    gauss[i]=fb->Integral(-lim+i,-lim+i+1);
    sum+=gauss[i]; 
  }
  for (int i=0; i<N; ++i)
    gauss[i] /= sum;

  for (int i=0; i<g->GetN(); ++i){
    double avy=0., avx=0., x, x0, y;
    int points=0;
    for (int j=i-N/2; j<=i+N/2; ++j){
      if      (j<0)          old->GetPoint(0, x, y);
      else if (j>=g->GetN()) old->GetPoint(old->GetN()-1, x, y);
      else                   old->GetPoint(j, x, y);
      if (i==j) x0=x;
      avy += y * gauss[points];
      avx += x * gauss[points];
      ++points;
    }
    if (i-N/2<0 || i+N/2>=g->GetN()) g->SetPoint(i, x0, avy);
    else g->SetPoint(i, avx, avy); 
  }
  delete old;
}


void Smooth2D(TGraph * g, int N)
{
  //TGraph * old = (TGraph*)g->Clone();
  TGraph * old = Close2D(g);
  //int N = (n%2==0?n+1:n);
  if (N>2*g->GetN()) N=2*g->GetN()-1;


  double gauss[N];
  double sigma = (double)N/4.;
  double sum=0;
  double lim=(double)N/2.;
  TF1 *fb = new TF1("fb","gaus(0)",-lim,lim);
  fb->SetParameter(0, 1./(sqrt(2*3.1415)*sigma) );
  fb->SetParameter(1, 0);
  fb->SetParameter(2, sigma);
  for (int i=0; i<N; ++i){
    gauss[i]=fb->Integral(-lim+i,-lim+i+1);
    sum+=gauss[i]; 
  }
  for (int i=0; i<N; ++i)
    gauss[i] /= sum;

  for (int i=0; i<g->GetN(); ++i){
    double avy=0., avx=0, x, x0, y;
    int points=0;
    for (int j=i-N/2; j<=i+N/2; ++j){
      //if      (j<0)          old->GetPoint(old->GetN()+j, x, y);
      //else if (j>=g->GetN()) old->GetPoint(j-old->GetN(), x, y);
      if      (j<0)          old->GetPoint(0, x, y);
      else if (j>=g->GetN()) old->GetPoint(old->GetN()-1, x, y);
      else                   old->GetPoint(j, x, y);
      if (i==j) x0=x;
      avy += y * gauss[points];
      avx += x * gauss[points];
      ++points;
    }
    g->SetPoint(i, avx, avy); 
  }
  delete old;
}

TGraph * Close2D(TGraph * g)
{
  double x, y;
  g->GetPoint(0,x,y);
  g->SetPoint(g->GetN(),x,y);
  
  int i=0;
  for (;i<g->GetN();++i){
    g->GetPoint(i,x,y);
    if (x<450&&y<450) break;
  }  
  int p=0;
  TGraph * f = new TGraph(0);
  for (int j=i+1;j!=i;++j){
//  for (int j=i+1;p<=g->GetN();++j){
    if (j>=g->GetN()) j=0;
    g->GetPoint(j,x,y);
    if (y<110+(x-120)*390/442||(x<330&&y<1000)||(x<500&&y<600)) continue;
    f->SetPoint(p++,x,y);
  }  
  return f; 
}


template class PlotTools<SusyScan>;
template class PlotTools<GeneratorMasses>;
Revision:	1.12
Committed:	Wed Jun 22 15:03:36 2011 UTC (13 years, 10 months ago) by auterman
Content type:	text/plain
Branch:	MAIN
CVS Tags:	JHEP2010, HEAD
Changes since 1.11:	+3 -0 lines
Log Message:	2010 RA2 paper
#	Content
1	#include "PlotTools.h"
2	#include "SusyScan.h"
3	#include "GeneratorMasses.h"
4
5	#include <cmath>
6	#include <algorithm>
7	#include <iostream>
8
9	#include "TGraph.h"
10	#include "TF1.h"
11	#include "TH2.h"
12	#include "TH2F.h"
13	#include "TObjArray.h"
14	#include "TPad.h"
15	#include "TCanvas.h"
16	#include "TRint.h"
17	#include "TROOT.h"
18
19	template<class T>
20	TGraph * PlotTools<T>::Line( double(x)(const T), double(y)(const T),
21	double(func)(const T), const double mass, const double diff )
22	{
23	TGraph * result = new TGraph(1);
24	std::sort(_scan->begin(),_scan->end(),sort_by(x));
25	int i=0;
26	for (typename std::vector<T*>::const_iterator it=_scan->begin();it!=_scan->end();++it){
27	if ( fabs(func( *it)-mass)<diff )
28	result->SetPoint(i++, x(it), y(it));
29	}
30	result->SetLineWidth(0.5);
31	result->SetLineColor(kGray);
32	return result;
33	}
34
35
36	template<class T>
37	void PlotTools<T>::Area( TH2h, double(x)(const T), double(y)(const T*),
38	double(f)(const T) )
39	{
40	//std::cout<< _scan->size() <<", &scan="<<_scan<<std::endl;
41	for (typename std::vector<T*>::const_iterator it=_scan->begin();it!=_scan->end();++it){
42	h->SetBinContent( h->GetXaxis()->FindBin(x(*it)),
43	h->GetYaxis()->FindBin(y(it)), f(it) );
44
45	// std::cout <<"m0="<< (it)->Mzero<<", m1/2="<< (it)->Mhalf<< ", sq="<<(it)->MUL << ", gl=" <<(it)->MGL <<std::endl;
46
47	}
48	}
49
50	template<class T>
51	void PlotTools<T>::Graph( TGraphg, double(x)(const T), double(y)(const T*),double ymin)
52	{
53	unsigned i = g->GetN();
54	std::sort(_scan->begin(),_scan->end(),sort_by(x));
55	for (typename std::vector<T*>::const_iterator it=_scan->begin();it!=_scan->end();++it){
56	if (y(it)>=ymin) g->SetPoint(i++, x(it), y(*it));
57	//std::cout << i << ": x="<<x(it)<< ", y="<<y(it)<< std::endl;
58	}
59	}
60
61	template<class T>
62	std::vector<TGraph> PlotTools<T>::GetContours005(TH2h, int ncont)
63	{
64	double conts[1];
65	conts[0] = 0.05;
66	if (!h) return std::vector<TGraph*>();
67	TH2 * plot = (TH2*)h->Clone();
68	plot->SetContour(1, conts);
69	plot->SetFillColor(1);
70	plot->Draw("CONT Z List");
71	gPad->Update();
72	TObjArray contours = (TObjArray)gROOT->GetListOfSpecials()->FindObject("contours");
73	int ncontours = contours->GetSize();
74	std::vector<TGraph*> result;
75	for (int i=0;i<ncontours;++i){
76	TList list = (TList)contours->At(i);
77	TGraph* curv = (TGraph*)list->First();
78	if (curv) result.push_back( curv );
79	for(int j = 0; j < list->GetSize(); j++){
80	curv = (TGraph*)list->After(curv); // Get Next graph
81	if (curv) result.push_back( curv );
82	}
83	}
84	delete plot;
85	std::sort(result.begin(),result.end(),sort_TGraph());
86	return result;
87	}
88
89	template<class T>
90	TGraph * PlotTools<T>::GetContour005(TH2*h, int ncont, int flag)
91	{
92	return (TGraph*)GetContours005(h, ncont).at(flag)->Clone();
93	}
94
95	template<class T>
96	std::vector<TGraph> PlotTools<T>::GetContours(TH2h, int ncont)
97	{
98	if (!h) return std::vector<TGraph*>();
99	TH2 * plot = (TH2*)h->Clone();
100	plot->SetContour(ncont);
101	plot->SetFillColor(1);
102	plot->Draw("CONT Z List");
103	gPad->Update();
104	TObjArray contours = (TObjArray)gROOT->GetListOfSpecials()->FindObject("contours");
105	int ncontours = contours->GetSize();
106	std::vector<TGraph*> result;
107	for (int i=0;i<ncontours;++i){
108	TList list = (TList)contours->At(i);
109	TGraph* curv = (TGraph*)list->First();
110	if (curv) result.push_back( curv );
111	for(int j = 0; j < list->GetSize(); j++){
112	curv = (TGraph*)list->After(curv); // Get Next graph
113	if (curv) result.push_back( curv );
114	}
115	}
116	delete plot;
117	std::sort(result.begin(),result.end(),sort_TGraph());
118	return result;
119	}
120
121	template<class T>
122	TGraph * PlotTools<T>::GetContour(TH2*h, int ncont, int flag)
123	{
124	return (TGraph*)GetContours(h, ncont).at(flag)->Clone();
125	}
126
127	template<class T>
128	TGraph * PlotTools<T>::GetContour(TH2h,double(x)(const T), double(y)(const T*),
129	double(func)(const T), int ncont, int flag,
130	int color, int style){
131	TH2hist = (TH2)h->Clone();
132	Area(hist, x, y, func);
133	TGraph * graph = GetContour(hist, ncont, flag);
134	graph->SetLineColor(color);
135	graph->SetLineStyle(style);
136	return graph;
137	}
138
139	template<class T>
140	TGraph * PlotTools<T>::GetContour005(TH2h,double(x)(const T), double(y)(const T*),
141	double(func)(const T), int ncont, int flag,
142	int color, int style){
143	TH2hist = (TH2)h->Clone();
144	Area(hist, x, y, func);
145	TGraph * graph = GetContour005(hist, ncont, flag);
146	graph->SetLineColor(color);
147	graph->SetLineStyle(style);
148	return graph;
149	}
150
151	template<class T>
152	void PlotTools<T>::Print(double(f)(const T), double(x)(const T), double(y)(const T), TGraph*g, double p)
153	{
154	for (typename std::vector<T*>::const_iterator it=_scan->begin();it!=_scan->end();++it){
155	for (int j=0; j<g->GetN(); ++j) {
156	double gx, gy;
157	g->GetPoint(j, gx, gy);
158	if ( (x(it)-gx)(x(it)-gx) + (y(it)-gy)(y(it)-gy) < p*p)
159	std::cout << x(it) << ", " << y(it)
160	<< " :: " << f(*it) << std::endl;
161	}
162	}
163
164	}
165	template<class T>
166	void PlotTools<T>::Print(double(f)(const T),
167	double(x1)(const T), double(x2)(const T),
168	double(x3)(const T), double(x4)(const T),
169	TGraph*g, double p)
170	{
171	for (typename std::vector<T*>::const_iterator it=_scan->begin();it!=_scan->end();++it){
172	for (int j=0; j<g->GetN(); ++j) {
173	double gx, gy;
174	g->GetPoint(j, gx, gy);
175	if ( (x1(it)-gx)(x1(it)-gx) + (x2(it)-gy)(x2(it)-gy) < p*p)
176	std::cout << x1(it) << ", " << x2(it) << ", "
177	<< x3(it) << ", " << x4(it)
178	<< " :: " << f(*it) << std::endl;
179	}
180	}
181
182	}
183
184	template<class T>
185	TGraph * PlotTools<T>::ChooseBest(TGraphr1,TGraphr2,TGraphg1,TGraphg2,double x,double y)
186	{
187	TGraph * res = new TGraph(0);
188	for (int i=0; i<r1->GetN(); ++i) {
189	double rx, r1y;
190	r1->GetPoint(i,rx,r1y);
191	double r2y = r2->Eval(rx);
192	double g1y=g1->Eval(rx);
193	double resy = (r1y>=r2y?g1y:g2->Eval(rx));
194	if (rx>x&&r1y>y)
195	res->SetPoint(i,rx,resy);
196	else
197	res->SetPoint(i,rx,g1y);
198	}
199	res->SetLineColor(g1->GetLineColor());
200	res->SetLineWidth(g1->GetLineWidth());
201	res->SetLineStyle(g1->GetLineStyle());
202	res->SetFillColor(g1->GetFillColor());
203	res->SetFillStyle(g1->GetFillStyle());
204	return res;
205	}
206
207	template<class T>
208	TH2 * PlotTools<T>::BinWiseOr(TH2h1, TH2h2)
209	{
210	TH2 * res = (TH2*)h1->Clone();
211	for (int x=0; x<=res->GetXaxis()->GetNbins(); ++x)
212	for (int y=0; y<=res->GetYaxis()->GetNbins(); ++y)
213	if (h2->GetBinContent(x,y)>0.5)
214	res->SetBinContent(x,y,h2->GetBinContent(x,y));
215	return res;
216	}
217
218	template<class T>
219	TGraph * PlotTools<T>::ModifyExpSigma(TGraphg1, TGraphg2, TGraph*g3)
220	{
221	TGraph * res = new TGraph(0);
222	for (int i=0; i<g1->GetN(); ++i) {
223	double x, y;
224	g1->GetPoint(i,x,y);
225	double y2 = g2->Eval(x);
226	double y3 = g3->Eval(x);
227	res->SetPoint(i,x,y-y2+y3);
228	}
229	res->SetLineColor(g1->GetLineColor());
230	res->SetLineStyle(g1->GetLineStyle());
231	res->SetLineWidth(g1->GetLineWidth());
232	res->SetFillColor(g1->GetFillColor());
233	res->SetFillStyle(g1->GetFillStyle());
234	return res;
235	}
236
237	double GetX(TGraph* g, double x0, double y0){
238	double dn=999999999, up=9999999999, minx, maxx;
239	for (int i=0; i<g->GetN(); ++i) {
240	double x, y;
241	g->GetPoint(i,x,y);
242	if (y==y0) return x;
243	if ( y<y0 && (x-x0)(x-x0)+(y-y0)(y-y0)<dn ){
244	dn=(x-x0)(x-x0)+(y-y0)(y-y0);
245	minx=x;
246	}
247	if ( y>y0 && (x-x0)(x-x0)+(y-y0)(y-y0)<up ){
248	up=(x-x0)(x-x0)+(y-y0)(y-y0);
249	maxx=x;
250	}
251	}
252	return (minx+maxx)/2.;
253	}
254
255	template<class T>
256	TGraph * PlotTools<T>::ModifyExpSigmaY(TGraphg1, TGraphg2, TGraph*g3)
257	{
258	TGraph * res = new TGraph(0);
259	for (int i=0; i<g1->GetN(); ++i) {
260	double x, y;
261	g1->GetPoint(i,x,y);
262	double x2,y2,x3,y3;
263	x2=GetX(g2,x,y);
264	x3=GetX(g3,x,y);
265	res->SetPoint(i,x-x2+x3,y);
266	}
267	res->SetLineColor(g1->GetLineColor());
268	res->SetLineStyle(g1->GetLineStyle());
269	res->SetLineWidth(g1->GetLineWidth());
270	res->SetFillColor(g1->GetFillColor());
271	res->SetFillStyle(g1->GetFillStyle());
272	return res;
273	}
274
275
276	template<class T>
277	bool PlotTools<T>::sort_TGraph::operator()(const TGraphg1, const TGraphg2)
278	{
279	return g1->GetN() > g2->GetN();
280	}
281
282	TGraph * MakeBand(TGraph g1, TGraph g2, bool b){
283	TGraph * res = new TGraph(g1->GetN()+g2->GetN());
284	int p=0;
285	for (int i=0; i<g1->GetN(); ++i) {
286	double x, y;
287	g1->GetPoint(i, x, y);
288	res->SetPoint(p++, x, y);
289	}
290	for (int i=g2->GetN()-1; i>=0; --i) {
291	double x, y;
292	g2->GetPoint(i, x, y);
293	res->SetPoint(p++, x, y);
294	}
295	res->SetLineColor( g1->GetLineColor() );
296	res->SetFillColor( g2->GetLineColor() );
297	res->SetFillStyle(4050);
298	return res;
299	}
300
301	void Smooth(TGraph * g, int N)
302	{
303	TGraph * old = (TGraph*)g->Clone();
304	//int N = (n%2==0?n+1:n);
305	if (N>2g->GetN()) N=2g->GetN()-1;
306
307
308	double gauss[N];
309	double sigma = (double)N/4.;
310	double sum=0;
311	double lim=(double)N/2.;
312	TF1 *fb = new TF1("fb","gaus(0)",-lim,lim);
313	fb->SetParameter(0, 1./(sqrt(23.1415)sigma) );
314	fb->SetParameter(1, 0);
315	fb->SetParameter(2, sigma);
316	for (int i=0; i<N; ++i){
317	gauss[i]=fb->Integral(-lim+i,-lim+i+1);
318	sum+=gauss[i];
319	}
320	for (int i=0; i<N; ++i)
321	gauss[i] /= sum;
322
323	for (int i=0; i<g->GetN(); ++i){
324	double avy=0., avx=0., x, x0, y;
325	int points=0;
326	for (int j=i-N/2; j<=i+N/2; ++j){
327	if (j<0) old->GetPoint(0, x, y);
328	else if (j>=g->GetN()) old->GetPoint(old->GetN()-1, x, y);
329	else old->GetPoint(j, x, y);
330	if (i==j) x0=x;
331	avy += y * gauss[points];
332	avx += x * gauss[points];
333	++points;
334	}
335	if (i-N/2<0 \|\| i+N/2>=g->GetN()) g->SetPoint(i, x0, avy);
336	else g->SetPoint(i, avx, avy);
337	}
338	delete old;
339	}
340
341
342	void Smooth2D(TGraph * g, int N)
343	{
344	//TGraph * old = (TGraph*)g->Clone();
345	TGraph * old = Close2D(g);
346	//int N = (n%2==0?n+1:n);
347	if (N>2g->GetN()) N=2g->GetN()-1;
348
349
350	double gauss[N];
351	double sigma = (double)N/4.;
352	double sum=0;
353	double lim=(double)N/2.;
354	TF1 *fb = new TF1("fb","gaus(0)",-lim,lim);
355	fb->SetParameter(0, 1./(sqrt(23.1415)sigma) );
356	fb->SetParameter(1, 0);
357	fb->SetParameter(2, sigma);
358	for (int i=0; i<N; ++i){
359	gauss[i]=fb->Integral(-lim+i,-lim+i+1);
360	sum+=gauss[i];
361	}
362	for (int i=0; i<N; ++i)
363	gauss[i] /= sum;
364
365	for (int i=0; i<g->GetN(); ++i){
366	double avy=0., avx=0, x, x0, y;
367	int points=0;
368	for (int j=i-N/2; j<=i+N/2; ++j){
369	//if (j<0) old->GetPoint(old->GetN()+j, x, y);
370	//else if (j>=g->GetN()) old->GetPoint(j-old->GetN(), x, y);
371	if (j<0) old->GetPoint(0, x, y);
372	else if (j>=g->GetN()) old->GetPoint(old->GetN()-1, x, y);
373	else old->GetPoint(j, x, y);
374	if (i==j) x0=x;
375	avy += y * gauss[points];
376	avx += x * gauss[points];
377	++points;
378	}
379	g->SetPoint(i, avx, avy);
380	}
381	delete old;
382	}
383
384	TGraph * Close2D(TGraph * g)
385	{
386	double x, y;
387	g->GetPoint(0,x,y);
388	g->SetPoint(g->GetN(),x,y);
389
390	int i=0;
391	for (;i<g->GetN();++i){
392	g->GetPoint(i,x,y);
393	if (x<450&&y<450) break;
394	}
395	int p=0;
396	TGraph * f = new TGraph(0);
397	for (int j=i+1;j!=i;++j){
398	// for (int j=i+1;p<=g->GetN();++j){
399	if (j>=g->GetN()) j=0;
400	g->GetPoint(j,x,y);
401	if (y<110+(x-120)*390/442\|\|(x<330&&y<1000)\|\|(x<500&&y<600)) continue;
402	f->SetPoint(p++,x,y);
403	}
404	return f;
405	}
406
407
408	template class PlotTools<SusyScan>;
409	template class PlotTools<GeneratorMasses>;