| 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( TH2*h, 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( TGraph*g, 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>::GetContours005(TH2*h, 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 |
std::vector<TGraph*> PlotTools<T>::GetContours(TH2*h, int ncont)
|
| 94 |
{
|
| 95 |
if (!h) return std::vector<TGraph*>();
|
| 96 |
TH2 * plot = (TH2*)h->Clone();
|
| 97 |
plot->SetContour(ncont);
|
| 98 |
plot->SetFillColor(1);
|
| 99 |
plot->Draw("CONT Z List");
|
| 100 |
gPad->Update();
|
| 101 |
TObjArray *contours = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
|
| 102 |
int ncontours = contours->GetSize();
|
| 103 |
std::vector<TGraph*> result;
|
| 104 |
for (int i=0;i<ncontours;++i){
|
| 105 |
TList *list = (TList*)contours->At(i);
|
| 106 |
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 |
}
|
| 113 |
delete plot;
|
| 114 |
std::sort(result.begin(),result.end(),sort_TGraph());
|
| 115 |
return result;
|
| 116 |
}
|
| 117 |
|
| 118 |
template<class T>
|
| 119 |
TGraph * PlotTools<T>::GetContour(TH2*h, int ncont, int flag)
|
| 120 |
{
|
| 121 |
return (TGraph*)GetContours(h, ncont).at(flag)->Clone();
|
| 122 |
}
|
| 123 |
|
| 124 |
template<class T>
|
| 125 |
TGraph * PlotTools<T>::GetContour(TH2*h,double(*x)(const T*), double(*y)(const T*),
|
| 126 |
double(*func)(const T*), int ncont, int flag,
|
| 127 |
int color, int style){
|
| 128 |
TH2*hist = (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 |
|
| 136 |
template<class T>
|
| 137 |
TGraph * PlotTools<T>::GetContour005(TH2*h,double(*x)(const T*), double(*y)(const T*),
|
| 138 |
double(*func)(const T*), int ncont, int flag,
|
| 139 |
int color, int style){
|
| 140 |
TH2*hist = (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 |
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 |
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 |
}
|
| 177 |
}
|
| 178 |
|
| 179 |
}
|
| 180 |
|
| 181 |
template<class T>
|
| 182 |
TGraph * PlotTools<T>::ChooseBest(TGraph*r1,TGraph*r2,TGraph*g1,TGraph*g2,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(TH2*h1, TH2*h2)
|
| 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 |
|
| 215 |
template<class T>
|
| 216 |
TGraph * PlotTools<T>::ModifyExpSigma(TGraph*g1, TGraph*g2, 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(TGraph*g1, TGraph*g2, 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 |
|
| 273 |
template<class T>
|
| 274 |
bool PlotTools<T>::sort_TGraph::operator()(const TGraph*g1, const TGraph*g2)
|
| 275 |
{
|
| 276 |
return g1->GetN() > g2->GetN();
|
| 277 |
}
|
| 278 |
|
| 279 |
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 |
res->SetFillStyle(4050);
|
| 295 |
return res;
|
| 296 |
}
|
| 297 |
|
| 298 |
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>2*g->GetN()) N=2*g->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(2*3.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 |
double avy=0., avx=0., x, x0, y;
|
| 322 |
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 |
avy += y * gauss[points];
|
| 329 |
avx += x * gauss[points];
|
| 330 |
++points;
|
| 331 |
}
|
| 332 |
if (i-N/2<0 || i+N/2>=g->GetN()) g->SetPoint(i, x0, avy);
|
| 333 |
else g->SetPoint(i, avx, avy);
|
| 334 |
}
|
| 335 |
delete old;
|
| 336 |
}
|
| 337 |
|
| 338 |
|
| 339 |
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>2*g->GetN()) N=2*g->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(2*3.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 |
// for (int j=i+1;p<=g->GetN();++j){
|
| 396 |
if (j>=g->GetN()) j=0;
|
| 397 |
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 |
template class PlotTools<SusyScan>;
|
| 406 |
template class PlotTools<GeneratorMasses>;
|
