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) );
  } 
}

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