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
#	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	auterman	1.12
45			// std::cout <<"m0="<< (it)->Mzero<<", m1/2="<< (it)->Mhalf<< ", sq="<<(it)->MUL << ", gl=" <<(it)->MGL <<std::endl;
46
47	auterman	1.1	}
48			}
49
50			template<class T>
51	auterman	1.5	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	auterman	1.10	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	auterman	1.3	std::vector<TGraph> PlotTools<T>::GetContours(TH2h, int ncont)
97	auterman	1.1	{
98	auterman	1.3	if (!h) return std::vector<TGraph*>();
99			TH2 * plot = (TH2*)h->Clone();
100			plot->SetContour(ncont);
101	auterman	1.1	plot->SetFillColor(1);
102	auterman	1.4	plot->Draw("CONT Z List");
103	auterman	1.1	gPad->Update();
104			TObjArray contours = (TObjArray)gROOT->GetListOfSpecials()->FindObject("contours");
105			int ncontours = contours->GetSize();
106	auterman	1.3	std::vector<TGraph*> result;
107			for (int i=0;i<ncontours;++i){
108			TList list = (TList)contours->At(i);
109	auterman	1.4	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	auterman	1.3	}
116			delete plot;
117	auterman	1.5	std::sort(result.begin(),result.end(),sort_TGraph());
118	auterman	1.3	return result;
119			}
120
121			template<class T>
122			TGraph * PlotTools<T>::GetContour(TH2*h, int ncont, int flag)
123			{
124	auterman	1.4	return (TGraph*)GetContours(h, ncont).at(flag)->Clone();
125			}
126
127			template<class T>
128	auterman	1.6	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	auterman	1.7
139			template<class T>
140	auterman	1.10	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	auterman	1.7	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	auterman	1.11	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	auterman	1.7	}
180			}
181
182			}
183
184	auterman	1.9	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	auterman	1.7
218	auterman	1.11	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	auterman	1.6
276			template<class T>
277	auterman	1.5	bool PlotTools<T>::sort_TGraph::operator()(const TGraphg1, const TGraphg2)
278			{
279			return g1->GetN() > g2->GetN();
280	auterman	1.1	}
281	auterman	1.5
282	auterman	1.6	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	auterman	1.9	res->SetFillStyle(4050);
298	auterman	1.6	return res;
299			}
300
301	auterman	1.8	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	auterman	1.10	double avy=0., avx=0., x, x0, y;
325	auterman	1.8	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	auterman	1.10	avy += y * gauss[points];
332			avx += x * gauss[points];
333	auterman	1.8	++points;
334			}
335	auterman	1.10	if (i-N/2<0 \|\| i+N/2>=g->GetN()) g->SetPoint(i, x0, avy);
336			else g->SetPoint(i, avx, avy);
337	auterman	1.8	}
338			delete old;
339			}
340
341	auterman	1.1
342	auterman	1.10	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	auterman	1.11	// for (int j=i+1;p<=g->GetN();++j){
399			if (j>=g->GetN()) j=0;
400	auterman	1.10	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	auterman	1.1	template class PlotTools<SusyScan>;
409			template class PlotTools<GeneratorMasses>;