SusyScan/PlotScript/scan.C

#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 "scan.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 = 200;
  min_x  = 0.0;
  max_y  = 2000.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_x,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_x,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_x,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;
}

MyTGraph * 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;
   
   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( "mass_scan.dat" );
   cout << "Read in " << points.size() << " mass points." <<endl;

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

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

   ///Where is the LSP charged?
   TH2F   * ch_LSP = Area<SUSY_POINT>( points, charged_LSP );
   MyTGraph * 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
   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();

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