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>::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>
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>
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>
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 av=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;
      av += y * gauss[points];
      ++points;
    }
    //g->SetPoint(i, x0, av/(double)points); 
    g->SetPoint(i, x0, av); 
  }
  delete old;
}


template class PlotTools<SusyScan>;
template class PlotTools<GeneratorMasses>;
Revision:	1.9
Committed:	Wed Mar 2 13:31:33 2011 UTC (14 years, 2 months ago) by auterman
Content type:	text/plain
Branch:	MAIN
Changes since 1.8:	+36 -0 lines
Log Message:	ARC review improvements
#	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	}
46
47	template<class T>
48	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	std::vector<TGraph> PlotTools<T>::GetContours(TH2h, int ncont)
60	{
61	if (!h) return std::vector<TGraph*>();
62	TH2 * plot = (TH2*)h->Clone();
63	plot->SetContour(ncont);
64	plot->SetFillColor(1);
65	plot->Draw("CONT Z List");
66	gPad->Update();
67	TObjArray contours = (TObjArray)gROOT->GetListOfSpecials()->FindObject("contours");
68	int ncontours = contours->GetSize();
69	std::vector<TGraph*> result;
70	for (int i=0;i<ncontours;++i){
71	TList list = (TList)contours->At(i);
72	TGraph* curv = (TGraph*)list->First();
73	if (curv) result.push_back( curv );
74	for(int j = 0; j < list->GetSize(); j++){
75	curv = (TGraph*)list->After(curv); // Get Next graph
76	if (curv) result.push_back( curv );
77	}
78	}
79	delete plot;
80	std::sort(result.begin(),result.end(),sort_TGraph());
81	return result;
82	}
83
84	template<class T>
85	TGraph * PlotTools<T>::GetContour(TH2*h, int ncont, int flag)
86	{
87	return (TGraph*)GetContours(h, ncont).at(flag)->Clone();
88	}
89
90	template<class T>
91	TGraph * PlotTools<T>::GetContour(TH2h,double(x)(const T), double(y)(const T*),
92	double(func)(const T), int ncont, int flag,
93	int color, int style){
94	TH2hist = (TH2)h->Clone();
95	Area(hist, x, y, func);
96	TGraph * graph = GetContour(hist, ncont, flag);
97	graph->SetLineColor(color);
98	graph->SetLineStyle(style);
99	return graph;
100	}
101
102	template<class T>
103	void PlotTools<T>::Print(double(f)(const T), double(x)(const T), double(y)(const T), TGraph*g, double p)
104	{
105	for (typename std::vector<T*>::const_iterator it=_scan->begin();it!=_scan->end();++it){
106	for (int j=0; j<g->GetN(); ++j) {
107	double gx, gy;
108	g->GetPoint(j, gx, gy);
109	if ( (x(it)-gx)(x(it)-gx) + (y(it)-gy)(y(it)-gy) < p*p)
110	std::cout << x(it) << ", " << y(it) << " :: " << f(*it) << std::endl;
111	}
112	}
113
114	}
115
116	template<class T>
117	TGraph * PlotTools<T>::ChooseBest(TGraphr1,TGraphr2,TGraphg1,TGraphg2,double x,double y)
118	{
119	TGraph * res = new TGraph(0);
120	for (int i=0; i<r1->GetN(); ++i) {
121	double rx, r1y;
122	r1->GetPoint(i,rx,r1y);
123	double r2y = r2->Eval(rx);
124	double g1y=g1->Eval(rx);
125	double resy = (r1y>=r2y?g1y:g2->Eval(rx));
126	if (rx>x&&r1y>y)
127	res->SetPoint(i,rx,resy);
128	else
129	res->SetPoint(i,rx,g1y);
130	}
131	res->SetLineColor(g1->GetLineColor());
132	res->SetLineWidth(g1->GetLineWidth());
133	res->SetLineStyle(g1->GetLineStyle());
134	res->SetFillColor(g1->GetFillColor());
135	res->SetFillStyle(g1->GetFillStyle());
136	return res;
137	}
138
139	template<class T>
140	TH2 * PlotTools<T>::BinWiseOr(TH2h1, TH2h2)
141	{
142	TH2 * res = (TH2*)h1->Clone();
143	for (int x=0; x<=res->GetXaxis()->GetNbins(); ++x)
144	for (int y=0; y<=res->GetYaxis()->GetNbins(); ++y)
145	if (h2->GetBinContent(x,y)>0.5)
146	res->SetBinContent(x,y,h2->GetBinContent(x,y));
147	return res;
148	}
149
150
151	template<class T>
152	bool PlotTools<T>::sort_TGraph::operator()(const TGraphg1, const TGraphg2)
153	{
154	return g1->GetN() > g2->GetN();
155	}
156
157	TGraph * MakeBand(TGraph g1, TGraph g2, bool b){
158	TGraph * res = new TGraph(g1->GetN()+g2->GetN());
159	int p=0;
160	for (int i=0; i<g1->GetN(); ++i) {
161	double x, y;
162	g1->GetPoint(i, x, y);
163	res->SetPoint(p++, x, y);
164	}
165	for (int i=g2->GetN()-1; i>=0; --i) {
166	double x, y;
167	g2->GetPoint(i, x, y);
168	res->SetPoint(p++, x, y);
169	}
170	res->SetLineColor( g1->GetLineColor() );
171	res->SetFillColor( g2->GetLineColor() );
172	res->SetFillStyle(4050);
173	return res;
174	}
175
176	void Smooth(TGraph * g, int N)
177	{
178	TGraph * old = (TGraph*)g->Clone();
179	//int N = (n%2==0?n+1:n);
180	if (N>2g->GetN()) N=2g->GetN()-1;
181
182
183	double gauss[N];
184	double sigma = (double)N/4.;
185	double sum=0;
186	double lim=(double)N/2.;
187	TF1 *fb = new TF1("fb","gaus(0)",-lim,lim);
188	fb->SetParameter(0, 1./(sqrt(23.1415)sigma) );
189	fb->SetParameter(1, 0);
190	fb->SetParameter(2, sigma);
191	for (int i=0; i<N; ++i){
192	gauss[i]=fb->Integral(-lim+i,-lim+i+1);
193	sum+=gauss[i];
194	}
195	for (int i=0; i<N; ++i)
196	gauss[i] /= sum;
197
198	for (int i=0; i<g->GetN(); ++i){
199	double av=0., x, x0, y;
200	int points=0;
201	for (int j=i-N/2; j<=i+N/2; ++j){
202	if (j<0) old->GetPoint(0, x, y);
203	else if (j>=g->GetN()) old->GetPoint(old->GetN()-1, x, y);
204	else old->GetPoint(j, x, y);
205	if (i==j) x0=x;
206	av += y * gauss[points];
207	++points;
208	}
209	//g->SetPoint(i, x0, av/(double)points);
210	g->SetPoint(i, x0, av);
211	}
212	delete old;
213	}
214
215
216	template class PlotTools<SusyScan>;
217	template class PlotTools<GeneratorMasses>;