SusyScan/PlotScript/scan.cc

#include <algorithm>
#include <map>
#include <stdlib.h>
#include <iostream>
#include "Riostream.h"
#include <cmath>

#include <TApplication.h>
#include <TPostScript.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <TArrow.h>
#include <TGaxis.h>
#include <TPad.h>
#include <TColor.h>
#include <TStyle.h>
#include <TLatex.h>
#include <TMarker.h>
#include "TRint.h"
#include "TROOT.h"
#include "TObjArray.h"

#include "scan.h"
#include "SusyScan.h"


//TROOT root("GUI", "GUI test environement");

using namespace std;


double total_xsec(const SUSY_XSECS point)
{ return point.total(); }


bool charged_LSP(const SUSY_POINT point)
{ return ( (fabs(point.MZ1) - point.MTAU1)>0.0 ); }

bool gluino(const SUSY_POINT point, const double mass)
{ return ( fabs(point.MGL - mass)<5.0 ); }

bool ul(const SUSY_POINT point, const double mass)
{ return ( fabs(point.MUL - mass)<5.0 ); }

bool neutralino_1(const SUSY_POINT point, const double mass)
{ return ( fabs(fabs(point.MZ1) - mass)<1.0 ); }

bool sq_equal_gluino(const SUSY_POINT point, const double treshold)
{ return ( fabs(point.MUL - point.MGL)<treshold ); }


bool xs_sq_equal_gluino(const SUSY_XSECS point, const double treshold)
{ if (point.sgsg()<1E-5&&point.sqsq()<1E-5) return false; else return (fabs(point.sgsg() - point.sqsq())<treshold*point.sqsq()/100. ); }


TScan::TScan()
{
  plot_id = 0;
  bins_x = 200;
  min_x  = 0.0;
  max_x  = 2000.0;
  bins_y = 40;
  min_y  = 0.0;
  max_y  = 1000.0;
}

TScan::~TScan()
{
}

template <class T> MyTGraph * TScan::IsoMassLine(std::vector<T> pnts, bool(*func)(T,double), double mass )
{
  MyTGraph * result = new MyTGraph(1);
  int i=0;
  for (typename std::vector<T>::const_iterator it=pnts.begin();
       it!=pnts.end(); ++it) {
     if ( func( *it, mass) )
       result->SetPoint(i++, it->MZERO, it->MHALF);
  } 
  return result;
}

template <class T> TH2F * TScan::Area(std::vector<T> pnts, double(*func)(T))
{
  char * name = new char[255];
  sprintf(name,"Area_plot%d",++plot_id);
  TH2F * result = new TH2F(name,"",bins_x,min_x,max_x,bins_y,min_y,max_y);
  for (typename std::vector<T>::const_iterator it=pnts.begin();
       it!=pnts.end(); ++it) {
     result->Fill(it->MZERO, it->MHALF, func(*it) );
  } 
  delete name;
  return result;
}

template <class T> TH2F * TScan::Area(std::vector<T> pnts, bool(*func)(T))
{
  char * name = new char[255];
  sprintf(name,"Area_plot%d",++plot_id);
  TH2F * result = new TH2F(name,"",bins_x,min_x,max_x,bins_y,min_y,max_y);
  for (typename std::vector<T>::const_iterator it=pnts.begin();
       it!=pnts.end(); ++it) {
     if ( func( *it ) ) 
       result->Fill(it->MZERO, it->MHALF, 1.);
  } 
  delete name;
  return result;
}

TH2F * TScan::Area51()
{
  char * name = new char[255];
  sprintf(name,"Area51_plot%d",++plot_id);
  TH2F * result = new TH2F(name,"",bins_x,min_x,max_x,bins_y,min_y,max_y);
  for (int x=1; x<=bins_x; ++x){
    for (int y=1; y<=bins_y; ++y){
      result->SetBinContent(x,y,1.);
    }  
  }
  for (std::vector<SUSY_POINT>::const_iterator it=points.begin();
       it!=points.end(); ++it) {
    result->SetBinContent((int)(bins_x*it->MZERO/(max_x-min_x)),
                          (int)(bins_y*it->MHALF/(max_y-min_y)),0.);
  }     
  delete name;
  return result;
}

TGraph * TScan::GetContour(TH2F*plot, int flag)
{
/*
   char * name = new char[255];
   sprintf(name,"canv%d",++plot_id);
   TCanvas * canv = new TCanvas(name,"c1",600,600);
   canv->cd();
   plot->SetContour(2);
   plot->SetContourLevel(0,0);
   plot->SetContourLevel(1,0.5);
   plot->SetFillColor(1);
   plot->Draw("CONT List");
   canv->Update();
   TObjArray *contours = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
   int ncontours = contours->GetSize();
   if (ncontours==0) return 0;
   //cout << "N contours = " << ncontours << endl;
   TList *list = (TList*)contours->At(0);
   MyTGraph *gr1 = (MyTGraph*)list->First();
   delete name;
   delete canv;
   delete list;
*/   
   plot->SetContour(2);
   plot->SetContourLevel(0,0);
   plot->SetContourLevel(1,0.5);
   plot->SetFillColor(1);
   plot->Draw("CONT List");
   gPad->Update();
   TObjArray *contours = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
   int ncontours      = contours->GetSize();
   cout << "N contours = " << ncontours << endl;
   TList *list        = (TList*)contours->At(0);
   //list->GetSize();
   TGraph *gr1 = (TGraph*)list->First();
   return gr1;
}

void TScan::ReadGeneratedMasses(std::string file)
{
   ifstream masses_file;
   masses_file.open(file.c_str());
   SUSY_POINT p;
   while (1) {
      masses_file >> p.MZERO      
                  >> p.MHALF      
                  >> p.TANB       
                  >> p.SGNMU      
                  >> p.AZERO      
                  >> p.MTOP       
                  >> p.muQ        
                  >> p.Q          
                  >> p.M1         
                  >> p.M2
                  >> p.M3         
                  >> p.MGL        
                  >> p.MUL        
                  >> p.MB1        
                  >> p.MSN        
                  >> p.MNTAU      
                  >> p.MZ1        
                  >> p.MW1        
                  >> p.MHL        
                  >> p.MUR        
                  >> p.MB2        
                  >> p.MEL        
                  >> p.MTAU1      
                  >> p.MZ2        
                  >> p.MW2        
                  >> p.MHH        
                  >> p.MDL        
                  >> p.MT1        
                  >> p.MER        
                  >> p.MTAU2      
                  >> p.MZ3        
                  >> p.MHA        
                  >> p.MDR        
                  >> p.MT2        
                  >> p.MZ4        
                  >> p.MHp;

      if (!masses_file.good()) break;
      if (fabs(p.SGNMU)!=1.) {
         cerr << "check lines near m0=" << p.MZERO << ", m1/2=" << p.MHALF << endl;
         break;
      } 
      points.push_back(p);
   }
}

void TScan::ReadGeneratedXsects(std::string file)
{
   ifstream xsects_file;
   xsects_file.open(file.c_str());
   SUSY_XSECS x;
   while (1) {
      xsects_file >> x.MZERO      
                  >> x.MHALF      
                  >> x.TANB       
                  >> x.SGNMU      
                  >> x.AZERO      
                  >> x.MTOP       
                  >> x.XS0
                  >> x.XS201
                  >> x.XS202
                  >> x.XS204
                  >> x.XS205
                  >> x.XS207
                  >> x.XS208
                  >> x.XS209
                  >> x.XS210
                  >> x.XS211
                  >> x.XS212
                  >> x.XS213
                  >> x.XS214
                  >> x.XS216
                  >> x.XS217
                  >> x.XS218
                  >> x.XS219
                  >> x.XS220
                  >> x.XS221
                  >> x.XS222
                  >> x.XS223
                  >> x.XS224
                  >> x.XS225
                  >> x.XS226
                  >> x.XS227
                  >> x.XS228
                  >> x.XS229
                  >> x.XS230
                  >> x.XS231
                  >> x.XS232
                  >> x.XS233
                  >> x.XS234
                  >> x.XS235
                  >> x.XS236
                  >> x.XS237
                  >> x.XS238
                  >> x.XS239
                  >> x.XS240
                  >> x.XS241
                  >> x.XS242
                  >> x.XS243
                  >> x.XS244
                  >> x.XS246
                  >> x.XS247
                  >> x.XS248
                  >> x.XS249
                  >> x.XS250
                  >> x.XS251
                  >> x.XS252
                  >> x.XS253
                  >> x.XS254
                  >> x.XS256
                  >> x.XS258
                  >> x.XS259
                  >> x.XS261
                  >> x.XS262
                  >> x.XS263
                  >> x.XS264
                  >> x.XS265
                  >> x.XS271
                  >> x.XS272
                  >> x.XS273
                  >> x.XS274
                  >> x.XS275
                  >> x.XS276
                  >> x.XS277
                  >> x.XS278
                  >> x.XS279
                  >> x.XS280
                  >> x.XS281
                  >> x.XS282
                  >> x.XS283
                  >> x.XS284
                  >> x.XS285
                  >> x.XS286
                  >> x.XS287
                  >> x.XS288
                  >> x.XS289
                  >> x.XS290
                  >> x.XS291
                  >> x.XS292
                  >> x.XS293
                  >> x.XS294
                  >> x.XS295
                  >> x.XS296
                  >> x.XS297
                  >> x.XS298
                  >> x.XS299
                  >> x.XS300
                  >> x.XS301;
      if (!xsects_file.good()) break;
      if (fabs(x.SGNMU)!=1.) {
         cerr << "check lines near m0=" << x.MZERO << ", m1/2=" << x.MHALF << endl;
         break;
      } 
      xsects.push_back(x);
   }
}

int TScan::DoStuff()
{
   ///Read in the generated masses
   ReadGeneratedMasses( "tb10_mu1_a0_massscan.dat" );
   cout << "Read in " << points.size() << " mass points." <<endl;

   ///Read in the generated x-sections
   ReadGeneratedXsects( "tb10mu1a0s7tev_scan_xsec.dat" );
   cout << "Read in " << xsects.size() << " x-section points." <<endl;

   ///Where are no mass-points generated?
   TH2F   * no_solution         = Area51(  );
   TGraph * contour_no_solution = GetContour( no_solution );
   if (contour_no_solution) contour_no_solution->SetLineColor(1);
   if (contour_no_solution) contour_no_solution->SetLineWidth(4);

   ///Where is the LSP charged?
   TH2F   * ch_LSP = Area<SUSY_POINT>( points, charged_LSP );
   TGraph * contour_charged_LSP = GetContour( ch_LSP );
   contour_charged_LSP->SetLineColor(6);
   contour_charged_LSP->SetLineWidth(4);

   ///Main Canvas
   TCanvas * c1 = new TCanvas("c1","c1",600,600);
   c1->SetFillStyle   ( 4000 );
   c1->SetLeftMargin  ( 0.15 );
   c1->SetRightMargin ( 0.15 );
   c1->SetBottomMargin( 0.10 );
   c1->cd();
   //c1->SetTopMargin   ( 0.05 );
   //TPostScript ps("plots.ps",-112);

   ///The base hist of the mSUGRA-plane with axis labels
   //TH2F * plane = new TH2F("plane","The mSUGRA parameter space;m_{0} [GeV];m_{1/2} [GeV]", bins_x,min_x,max_x,bins_y,min_x,max_y);
   SetBins(80,80);//xsec-binning is 25 GeV
   TH2F * plane = Area<SUSY_XSECS>( xsects, total_xsec );
   SetBins(400,400);//mass-binning is 5 GeV
   plane->SetTitle("The mSUGRA parameter space;m_{0} [GeV];m_{1/2} [GeV];#sigma_{SUSY incl.} [pb]");
   plane->GetYaxis()->SetTitleOffset(  1.8 );
   plane->GetZaxis()->SetTitleOffset(  1.1 );
   plane->GetZaxis()->SetRangeUser(0.00001, 1000.);
   gPad->SetLogz(1);
   plane->Draw("COLZ");

   ///Now, add what ever should be displayed else
   //---------------------------------------------------------------------------//
   TLatex t;
   t.SetTextColor(14);
   t.SetTextAlign(12);
   t.SetTextSize(0.02);
   
   ///gluino
   gStyle->SetLineColor(15);//gray
   MyTGraph * m_gluino5  = IsoMassLine( points, gluino,  500. ); m_gluino5->Draw();
   MyTGraph * m_gluino10 = IsoMassLine( points, gluino, 1000. ); m_gluino10->Draw();
   MyTGraph * m_gluino20 = IsoMassLine( points, gluino, 2000. ); m_gluino20->Draw();
   MyTGraph * m_gluino30 = IsoMassLine( points, gluino, 3000. ); m_gluino30->Draw();
   MyTGraph * m_gluino40 = IsoMassLine( points, gluino, 4000. ); m_gluino40->Draw();
   t.DrawLatex(360.,450.,"#tilde{g}=1000 GeV");
   t.DrawLatex(600.,930.,"#tilde{g}=2000 GeV");
   t.DrawLatex(870.,1420.,"#tilde{g}=3000 GeV");
   t.DrawLatex(1150.,1930.,"#tilde{g}=4000 GeV");
   
   ///neutralino
   //MyTGraph * m_neutralino1 = IsoMassLine( neutralino_1, 500. ); m_neutralino1->Draw();
   
   ///left-handed up-like squarks
   MyTGraph * m_ul5  = IsoMassLine( points, ul,  500. ); m_ul5->Draw();
   MyTGraph * m_ul10 = IsoMassLine( points, ul, 1000. ); m_ul10->Draw();
   MyTGraph * m_ul15 = IsoMassLine( points, ul, 1500. ); m_ul15->Draw();
   MyTGraph * m_ul20 = IsoMassLine( points, ul, 2000. ); m_ul20->Draw();
   MyTGraph * m_ul25 = IsoMassLine( points, ul, 2500. ); m_ul25->Draw();
   MyTGraph * m_ul30 = IsoMassLine( points, ul, 3000. ); m_ul30->Draw();
   MyTGraph * m_ul35 = IsoMassLine( points, ul, 3500. ); m_ul35->Draw();
   MyTGraph * m_ul40 = IsoMassLine( points, ul, 4000. ); m_ul40->Draw();
   t.DrawLatex(800.,250.,"#tilde{u}_{L}=1000 GeV");
   t.DrawLatex(1670.,470.,"#tilde{u}_{L}=2000 GeV");
   t.DrawLatex(1670.,1250.,"#tilde{u}_{L}=3000 GeV");
   t.DrawLatex(1670.,1930.,"#tilde{u}_{L}=4000 GeV");

   ///Contours
/*
   if (contour_no_solution) {
   //contour_no_solution->Add(max_x+5,430);
   //contour_no_solution->Add(max_x+5,min_y-5);
   //contour_no_solution->Add(min_x,min_y-5);
   contour_no_solution->SetFillColor( 1 );
   contour_no_solution->SetFillStyle( 3454 );
   contour_no_solution->Draw("lf");
   contour_no_solution->Draw();
   }
*/
   if (contour_charged_LSP){
   //contour_charged_LSP->Add(min_x-5, 726.25);
   //contour_charged_LSP->Add(min_x-5, max_y+5);
   contour_charged_LSP->SetFillColor( 6 );
   contour_charged_LSP->SetFillStyle( 3454 );
   contour_charged_LSP->SetLineColor( 6 );
   contour_charged_LSP->SetLineWidth( 4 );
   contour_charged_LSP->Draw("lf");
   contour_charged_LSP->Draw();
   }
   t.SetTextSize(0.03);t.SetTextColor( 6 );
   t.DrawLatex(200.,1800.,"charged LSP (#tilde{#tau})");

   // m(gluino)==m(squark)
   gStyle->SetLineColor(1);//black
   MyTGraph * m_ulgl  = IsoMassLine( points, sq_equal_gluino,  0.2 ); 
   m_ulgl->SetLineWidth(3); m_ulgl->Draw();
   t.SetTextSize(0.03);t.SetTextColor( 1 );
   t.DrawLatex( 750.,1200.,"m(#tilde{u}_{L}) < m(#tilde{g})");
   t.DrawLatex(1350.,1150.,"m(#tilde{u}_{L}) > m(#tilde{g})");


   // xs(gluino)==xs(squark)
   MyTGraph * xs_ulgl  = IsoMassLine( xsects, xs_sq_equal_gluino,  5.0 ); 
   xs_ulgl->SetLineWidth(3); xs_ulgl->Draw();
   t.SetTextSize(0.02);t.SetTextColor( 1 );
   t.DrawLatex( 890.,510.,"#sigma(#tilde{q}#tilde{q}) > #sigmam(#tilde{g}#tilde{g})");
   t.DrawLatex(1280.,450.,"#sigma(#tilde{q}#tilde{q}) < #sigmam(#tilde{g}#tilde{g})");

   TMarker m;
   m.SetMarkerStyle(29);
   m.SetMarkerColor(2);
   m.SetMarkerSize(1.5);
   //m.DrawMarker(1000,1000);  
   //DONE
   c1->SaveAs("dirt.eps");
   gApplication->Run();    //interactive

   //ps.NewPage();

}

int scan(void)
{
   gStyle->SetHistFillColor(0);
   gStyle->SetPalette(1);
   //gStyle->SetFillColor(0);
   //gStyle->SetOptStat(kFALSE);
   gStyle->SetCanvasColor(0);
   gStyle->SetCanvasBorderMode(0);
   gStyle->SetPadColor(0);
   gStyle->SetPadBorderMode(0);
   gStyle->SetFrameBorderMode(0);

   gStyle->SetTitleFillColor(0);
   gStyle->SetTitleBorderSize(0);
   gStyle->SetTitleX(0.10);
   gStyle->SetTitleY(0.98);
   gStyle->SetTitleW(0.8);
   gStyle->SetTitleH(0.06);

   gStyle->SetErrorX(0);
   gStyle->SetStatColor(0);
   gStyle->SetStatBorderSize(0);
   gStyle->SetStatX(0);
   gStyle->SetStatY(0);
   gStyle->SetStatW(0);
   gStyle->SetStatH(0);

   gStyle->SetTitleFont(22);
   gStyle->SetLabelFont(22,"X");
   gStyle->SetLabelFont(22,"Y");
   gStyle->SetLabelFont(22,"Z");
   gStyle->SetLabelSize(0.03,"X");
   gStyle->SetLabelSize(0.03,"Y");
   gStyle->SetLabelSize(0.03,"Z");
 
   //gROOT->SetStyle("MyStyle");

   
   TScan * scan = new TScan();
   scan->DoStuff();
   delete scan;
}

int main(int argc, char** argv)
{
  TRint *theApp = new TRint("ROOT example", &argc, argv);
  //theApp->Run();
  //TApplication theApp("App", &argc, argv);
  
  return scan();
}
Revision:	1.2
Committed:	Sat Nov 27 12:03:26 2010 UTC (14 years, 5 months ago) by auterman
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.1:	+0 -0 lines
State:	*FILE REMOVED*
Error occurred while calculating annotation data.
Log Message:	reorganization of code
#	Content
1	#include <algorithm>
2	#include <map>
3	#include <stdlib.h>
4	#include <iostream>
5	#include "Riostream.h"
6	#include <cmath>
7
8	#include <TApplication.h>
9	#include <TPostScript.h>
10	#include <TCanvas.h>
11	#include <TLegend.h>
12	#include <TArrow.h>
13	#include <TGaxis.h>
14	#include <TPad.h>
15	#include <TColor.h>
16	#include <TStyle.h>
17	#include <TLatex.h>
18	#include <TMarker.h>
19	#include "TRint.h"
20	#include "TROOT.h"
21	#include "TObjArray.h"
22
23	#include "scan.h"
24	#include "SusyScan.h"
25
26
27	//TROOT root("GUI", "GUI test environement");
28
29	using namespace std;
30
31
32	double total_xsec(const SUSY_XSECS point)
33	{ return point.total(); }
34
35
36	bool charged_LSP(const SUSY_POINT point)
37	{ return ( (fabs(point.MZ1) - point.MTAU1)>0.0 ); }
38
39	bool gluino(const SUSY_POINT point, const double mass)
40	{ return ( fabs(point.MGL - mass)<5.0 ); }
41
42	bool ul(const SUSY_POINT point, const double mass)
43	{ return ( fabs(point.MUL - mass)<5.0 ); }
44
45	bool neutralino_1(const SUSY_POINT point, const double mass)
46	{ return ( fabs(fabs(point.MZ1) - mass)<1.0 ); }
47
48	bool sq_equal_gluino(const SUSY_POINT point, const double treshold)
49	{ return ( fabs(point.MUL - point.MGL)<treshold ); }
50
51
52	bool xs_sq_equal_gluino(const SUSY_XSECS point, const double treshold)
53	{ if (point.sgsg()<1E-5&&point.sqsq()<1E-5) return false; else return (fabs(point.sgsg() - point.sqsq())<treshold*point.sqsq()/100. ); }
54
55
56	TScan::TScan()
57	{
58	plot_id = 0;
59	bins_x = 200;
60	min_x = 0.0;
61	max_x = 2000.0;
62	bins_y = 40;
63	min_y = 0.0;
64	max_y = 1000.0;
65	}
66
67	TScan::~TScan()
68	{
69	}
70
71	template <class T> MyTGraph * TScan::IsoMassLine(std::vector<T> pnts, bool(*func)(T,double), double mass )
72	{
73	MyTGraph * result = new MyTGraph(1);
74	int i=0;
75	for (typename std::vector<T>::const_iterator it=pnts.begin();
76	it!=pnts.end(); ++it) {
77	if ( func( *it, mass) )
78	result->SetPoint(i++, it->MZERO, it->MHALF);
79	}
80	return result;
81	}
82
83	template <class T> TH2F * TScan::Area(std::vector<T> pnts, double(*func)(T))
84	{
85	char * name = new char[255];
86	sprintf(name,"Area_plot%d",++plot_id);
87	TH2F * result = new TH2F(name,"",bins_x,min_x,max_x,bins_y,min_y,max_y);
88	for (typename std::vector<T>::const_iterator it=pnts.begin();
89	it!=pnts.end(); ++it) {
90	result->Fill(it->MZERO, it->MHALF, func(*it) );
91	}
92	delete name;
93	return result;
94	}
95
96	template <class T> TH2F * TScan::Area(std::vector<T> pnts, bool(*func)(T))
97	{
98	char * name = new char[255];
99	sprintf(name,"Area_plot%d",++plot_id);
100	TH2F * result = new TH2F(name,"",bins_x,min_x,max_x,bins_y,min_y,max_y);
101	for (typename std::vector<T>::const_iterator it=pnts.begin();
102	it!=pnts.end(); ++it) {
103	if ( func( *it ) )
104	result->Fill(it->MZERO, it->MHALF, 1.);
105	}
106	delete name;
107	return result;
108	}
109
110	TH2F * TScan::Area51()
111	{
112	char * name = new char[255];
113	sprintf(name,"Area51_plot%d",++plot_id);
114	TH2F * result = new TH2F(name,"",bins_x,min_x,max_x,bins_y,min_y,max_y);
115	for (int x=1; x<=bins_x; ++x){
116	for (int y=1; y<=bins_y; ++y){
117	result->SetBinContent(x,y,1.);
118	}
119	}
120	for (std::vector<SUSY_POINT>::const_iterator it=points.begin();
121	it!=points.end(); ++it) {
122	result->SetBinContent((int)(bins_x*it->MZERO/(max_x-min_x)),
123	(int)(bins_y*it->MHALF/(max_y-min_y)),0.);
124	}
125	delete name;
126	return result;
127	}
128
129	TGraph * TScan::GetContour(TH2F*plot, int flag)
130	{
131	/*
132	char * name = new char[255];
133	sprintf(name,"canv%d",++plot_id);
134	TCanvas * canv = new TCanvas(name,"c1",600,600);
135	canv->cd();
136	plot->SetContour(2);
137	plot->SetContourLevel(0,0);
138	plot->SetContourLevel(1,0.5);
139	plot->SetFillColor(1);
140	plot->Draw("CONT List");
141	canv->Update();
142	TObjArray contours = (TObjArray)gROOT->GetListOfSpecials()->FindObject("contours");
143	int ncontours = contours->GetSize();
144	if (ncontours==0) return 0;
145	//cout << "N contours = " << ncontours << endl;
146	TList list = (TList)contours->At(0);
147	MyTGraph gr1 = (MyTGraph)list->First();
148	delete name;
149	delete canv;
150	delete list;
151	*/
152	plot->SetContour(2);
153	plot->SetContourLevel(0,0);
154	plot->SetContourLevel(1,0.5);
155	plot->SetFillColor(1);
156	plot->Draw("CONT List");
157	gPad->Update();
158	TObjArray contours = (TObjArray)gROOT->GetListOfSpecials()->FindObject("contours");
159	int ncontours = contours->GetSize();
160	cout << "N contours = " << ncontours << endl;
161	TList list = (TList)contours->At(0);
162	//list->GetSize();
163	TGraph gr1 = (TGraph)list->First();
164	return gr1;
165	}
166
167	void TScan::ReadGeneratedMasses(std::string file)
168	{
169	ifstream masses_file;
170	masses_file.open(file.c_str());
171	SUSY_POINT p;
172	while (1) {
173	masses_file >> p.MZERO
174	>> p.MHALF
175	>> p.TANB
176	>> p.SGNMU
177	>> p.AZERO
178	>> p.MTOP
179	>> p.muQ
180	>> p.Q
181	>> p.M1
182	>> p.M2
183	>> p.M3
184	>> p.MGL
185	>> p.MUL
186	>> p.MB1
187	>> p.MSN
188	>> p.MNTAU
189	>> p.MZ1
190	>> p.MW1
191	>> p.MHL
192	>> p.MUR
193	>> p.MB2
194	>> p.MEL
195	>> p.MTAU1
196	>> p.MZ2
197	>> p.MW2
198	>> p.MHH
199	>> p.MDL
200	>> p.MT1
201	>> p.MER
202	>> p.MTAU2
203	>> p.MZ3
204	>> p.MHA
205	>> p.MDR
206	>> p.MT2
207	>> p.MZ4
208	>> p.MHp;
209
210	if (!masses_file.good()) break;
211	if (fabs(p.SGNMU)!=1.) {
212	cerr << "check lines near m0=" << p.MZERO << ", m1/2=" << p.MHALF << endl;
213	break;
214	}
215	points.push_back(p);
216	}
217	}
218
219	void TScan::ReadGeneratedXsects(std::string file)
220	{
221	ifstream xsects_file;
222	xsects_file.open(file.c_str());
223	SUSY_XSECS x;
224	while (1) {
225	xsects_file >> x.MZERO
226	>> x.MHALF
227	>> x.TANB
228	>> x.SGNMU
229	>> x.AZERO
230	>> x.MTOP
231	>> x.XS0
232	>> x.XS201
233	>> x.XS202
234	>> x.XS204
235	>> x.XS205
236	>> x.XS207
237	>> x.XS208
238	>> x.XS209
239	>> x.XS210
240	>> x.XS211
241	>> x.XS212
242	>> x.XS213
243	>> x.XS214
244	>> x.XS216
245	>> x.XS217
246	>> x.XS218
247	>> x.XS219
248	>> x.XS220
249	>> x.XS221
250	>> x.XS222
251	>> x.XS223
252	>> x.XS224
253	>> x.XS225
254	>> x.XS226
255	>> x.XS227
256	>> x.XS228
257	>> x.XS229
258	>> x.XS230
259	>> x.XS231
260	>> x.XS232
261	>> x.XS233
262	>> x.XS234
263	>> x.XS235
264	>> x.XS236
265	>> x.XS237
266	>> x.XS238
267	>> x.XS239
268	>> x.XS240
269	>> x.XS241
270	>> x.XS242
271	>> x.XS243
272	>> x.XS244
273	>> x.XS246
274	>> x.XS247
275	>> x.XS248
276	>> x.XS249
277	>> x.XS250
278	>> x.XS251
279	>> x.XS252
280	>> x.XS253
281	>> x.XS254
282	>> x.XS256
283	>> x.XS258
284	>> x.XS259
285	>> x.XS261
286	>> x.XS262
287	>> x.XS263
288	>> x.XS264
289	>> x.XS265
290	>> x.XS271
291	>> x.XS272
292	>> x.XS273
293	>> x.XS274
294	>> x.XS275
295	>> x.XS276
296	>> x.XS277
297	>> x.XS278
298	>> x.XS279
299	>> x.XS280
300	>> x.XS281
301	>> x.XS282
302	>> x.XS283
303	>> x.XS284
304	>> x.XS285
305	>> x.XS286
306	>> x.XS287
307	>> x.XS288
308	>> x.XS289
309	>> x.XS290
310	>> x.XS291
311	>> x.XS292
312	>> x.XS293
313	>> x.XS294
314	>> x.XS295
315	>> x.XS296
316	>> x.XS297
317	>> x.XS298
318	>> x.XS299
319	>> x.XS300
320	>> x.XS301;
321	if (!xsects_file.good()) break;
322	if (fabs(x.SGNMU)!=1.) {
323	cerr << "check lines near m0=" << x.MZERO << ", m1/2=" << x.MHALF << endl;
324	break;
325	}
326	xsects.push_back(x);
327	}
328	}
329
330	int TScan::DoStuff()
331	{
332	///Read in the generated masses
333	ReadGeneratedMasses( "tb10_mu1_a0_massscan.dat" );
334	cout << "Read in " << points.size() << " mass points." <<endl;
335
336	///Read in the generated x-sections
337	ReadGeneratedXsects( "tb10mu1a0s7tev_scan_xsec.dat" );
338	cout << "Read in " << xsects.size() << " x-section points." <<endl;
339
340	///Where are no mass-points generated?
341	TH2F * no_solution = Area51( );
342	TGraph * contour_no_solution = GetContour( no_solution );
343	if (contour_no_solution) contour_no_solution->SetLineColor(1);
344	if (contour_no_solution) contour_no_solution->SetLineWidth(4);
345
346	///Where is the LSP charged?
347	TH2F * ch_LSP = Area<SUSY_POINT>( points, charged_LSP );
348	TGraph * contour_charged_LSP = GetContour( ch_LSP );
349	contour_charged_LSP->SetLineColor(6);
350	contour_charged_LSP->SetLineWidth(4);
351
352	///Main Canvas
353	TCanvas * c1 = new TCanvas("c1","c1",600,600);
354	c1->SetFillStyle ( 4000 );
355	c1->SetLeftMargin ( 0.15 );
356	c1->SetRightMargin ( 0.15 );
357	c1->SetBottomMargin( 0.10 );
358	c1->cd();
359	//c1->SetTopMargin ( 0.05 );
360	//TPostScript ps("plots.ps",-112);
361
362	///The base hist of the mSUGRA-plane with axis labels
363	//TH2F * plane = new TH2F("plane","The mSUGRA parameter space;m_{0} [GeV];m_{1/2} [GeV]", bins_x,min_x,max_x,bins_y,min_x,max_y);
364	SetBins(80,80);//xsec-binning is 25 GeV
365	TH2F * plane = Area<SUSY_XSECS>( xsects, total_xsec );
366	SetBins(400,400);//mass-binning is 5 GeV
367	plane->SetTitle("The mSUGRA parameter space;m_{0} [GeV];m_{1/2} [GeV];#sigma_{SUSY incl.} [pb]");
368	plane->GetYaxis()->SetTitleOffset( 1.8 );
369	plane->GetZaxis()->SetTitleOffset( 1.1 );
370	plane->GetZaxis()->SetRangeUser(0.00001, 1000.);
371	gPad->SetLogz(1);
372	plane->Draw("COLZ");
373
374	///Now, add what ever should be displayed else
375	//---------------------------------------------------------------------------//
376	TLatex t;
377	t.SetTextColor(14);
378	t.SetTextAlign(12);
379	t.SetTextSize(0.02);
380
381	///gluino
382	gStyle->SetLineColor(15);//gray
383	MyTGraph * m_gluino5 = IsoMassLine( points, gluino, 500. ); m_gluino5->Draw();
384	MyTGraph * m_gluino10 = IsoMassLine( points, gluino, 1000. ); m_gluino10->Draw();
385	MyTGraph * m_gluino20 = IsoMassLine( points, gluino, 2000. ); m_gluino20->Draw();
386	MyTGraph * m_gluino30 = IsoMassLine( points, gluino, 3000. ); m_gluino30->Draw();
387	MyTGraph * m_gluino40 = IsoMassLine( points, gluino, 4000. ); m_gluino40->Draw();
388	t.DrawLatex(360.,450.,"#tilde{g}=1000 GeV");
389	t.DrawLatex(600.,930.,"#tilde{g}=2000 GeV");
390	t.DrawLatex(870.,1420.,"#tilde{g}=3000 GeV");
391	t.DrawLatex(1150.,1930.,"#tilde{g}=4000 GeV");
392
393	///neutralino
394	//MyTGraph * m_neutralino1 = IsoMassLine( neutralino_1, 500. ); m_neutralino1->Draw();
395
396	///left-handed up-like squarks
397	MyTGraph * m_ul5 = IsoMassLine( points, ul, 500. ); m_ul5->Draw();
398	MyTGraph * m_ul10 = IsoMassLine( points, ul, 1000. ); m_ul10->Draw();
399	MyTGraph * m_ul15 = IsoMassLine( points, ul, 1500. ); m_ul15->Draw();
400	MyTGraph * m_ul20 = IsoMassLine( points, ul, 2000. ); m_ul20->Draw();
401	MyTGraph * m_ul25 = IsoMassLine( points, ul, 2500. ); m_ul25->Draw();
402	MyTGraph * m_ul30 = IsoMassLine( points, ul, 3000. ); m_ul30->Draw();
403	MyTGraph * m_ul35 = IsoMassLine( points, ul, 3500. ); m_ul35->Draw();
404	MyTGraph * m_ul40 = IsoMassLine( points, ul, 4000. ); m_ul40->Draw();
405	t.DrawLatex(800.,250.,"#tilde{u}_{L}=1000 GeV");
406	t.DrawLatex(1670.,470.,"#tilde{u}_{L}=2000 GeV");
407	t.DrawLatex(1670.,1250.,"#tilde{u}_{L}=3000 GeV");
408	t.DrawLatex(1670.,1930.,"#tilde{u}_{L}=4000 GeV");
409
410	///Contours
411	/*
412	if (contour_no_solution) {
413	//contour_no_solution->Add(max_x+5,430);
414	//contour_no_solution->Add(max_x+5,min_y-5);
415	//contour_no_solution->Add(min_x,min_y-5);
416	contour_no_solution->SetFillColor( 1 );
417	contour_no_solution->SetFillStyle( 3454 );
418	contour_no_solution->Draw("lf");
419	contour_no_solution->Draw();
420	}
421	*/
422	if (contour_charged_LSP){
423	//contour_charged_LSP->Add(min_x-5, 726.25);
424	//contour_charged_LSP->Add(min_x-5, max_y+5);
425	contour_charged_LSP->SetFillColor( 6 );
426	contour_charged_LSP->SetFillStyle( 3454 );
427	contour_charged_LSP->SetLineColor( 6 );
428	contour_charged_LSP->SetLineWidth( 4 );
429	contour_charged_LSP->Draw("lf");
430	contour_charged_LSP->Draw();
431	}
432	t.SetTextSize(0.03);t.SetTextColor( 6 );
433	t.DrawLatex(200.,1800.,"charged LSP (#tilde{#tau})");
434
435	// m(gluino)==m(squark)
436	gStyle->SetLineColor(1);//black
437	MyTGraph * m_ulgl = IsoMassLine( points, sq_equal_gluino, 0.2 );
438	m_ulgl->SetLineWidth(3); m_ulgl->Draw();
439	t.SetTextSize(0.03);t.SetTextColor( 1 );
440	t.DrawLatex( 750.,1200.,"m(#tilde{u}_{L}) < m(#tilde{g})");
441	t.DrawLatex(1350.,1150.,"m(#tilde{u}_{L}) > m(#tilde{g})");
442
443
444	// xs(gluino)==xs(squark)
445	MyTGraph * xs_ulgl = IsoMassLine( xsects, xs_sq_equal_gluino, 5.0 );
446	xs_ulgl->SetLineWidth(3); xs_ulgl->Draw();
447	t.SetTextSize(0.02);t.SetTextColor( 1 );
448	t.DrawLatex( 890.,510.,"#sigma(#tilde{q}#tilde{q}) > #sigmam(#tilde{g}#tilde{g})");
449	t.DrawLatex(1280.,450.,"#sigma(#tilde{q}#tilde{q}) < #sigmam(#tilde{g}#tilde{g})");
450
451	TMarker m;
452	m.SetMarkerStyle(29);
453	m.SetMarkerColor(2);
454	m.SetMarkerSize(1.5);
455	//m.DrawMarker(1000,1000);
456	//DONE
457	c1->SaveAs("dirt.eps");
458	gApplication->Run(); //interactive
459
460	//ps.NewPage();
461
462	}
463
464	int scan(void)
465	{
466	gStyle->SetHistFillColor(0);
467	gStyle->SetPalette(1);
468	//gStyle->SetFillColor(0);
469	//gStyle->SetOptStat(kFALSE);
470	gStyle->SetCanvasColor(0);
471	gStyle->SetCanvasBorderMode(0);
472	gStyle->SetPadColor(0);
473	gStyle->SetPadBorderMode(0);
474	gStyle->SetFrameBorderMode(0);
475
476	gStyle->SetTitleFillColor(0);
477	gStyle->SetTitleBorderSize(0);
478	gStyle->SetTitleX(0.10);
479	gStyle->SetTitleY(0.98);
480	gStyle->SetTitleW(0.8);
481	gStyle->SetTitleH(0.06);
482
483	gStyle->SetErrorX(0);
484	gStyle->SetStatColor(0);
485	gStyle->SetStatBorderSize(0);
486	gStyle->SetStatX(0);
487	gStyle->SetStatY(0);
488	gStyle->SetStatW(0);
489	gStyle->SetStatH(0);
490
491	gStyle->SetTitleFont(22);
492	gStyle->SetLabelFont(22,"X");
493	gStyle->SetLabelFont(22,"Y");
494	gStyle->SetLabelFont(22,"Z");
495	gStyle->SetLabelSize(0.03,"X");
496	gStyle->SetLabelSize(0.03,"Y");
497	gStyle->SetLabelSize(0.03,"Z");
498
499	//gROOT->SetStyle("MyStyle");
500
501
502	TScan * scan = new TScan();
503	scan->DoStuff();
504	delete scan;
505	}
506
507	int main(int argc, char** argv)
508	{
509	TRint *theApp = new TRint("ROOT example", &argc, argv);
510	//theApp->Run();
511	//TApplication theApp("App", &argc, argv);
512
513	return scan();
514	}