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
#	User	Rev	Content
1	auterman	1.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	auterman	1.1.1.2	#include <TLatex.h>
18			#include <TMarker.h>
19	auterman	1.1
20
21			#include "scan.h"
22
23
24			TROOT root("GUI", "GUI test environement");
25
26			using namespace std;
27
28	auterman	1.1.1.2
29			double total_xsec(const SUSY_XSECS point)
30			{ return point.total(); }
31
32
33	auterman	1.1	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	auterman	1.1.1.2	{ return ( fabs(point.MGL - mass)<5.0 ); }
38	auterman	1.1
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	auterman	1.1.1.2	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	auterman	1.1
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	auterman	1.1.1.2	template <class T> MyTGraph * TScan::IsoMassLine(std::vector<T> pnts, bool(*func)(T,double), double mass )
69	auterman	1.1	{
70			MyTGraph * result = new MyTGraph(1);
71			int i=0;
72	auterman	1.1.1.2	for (typename std::vector<T>::const_iterator it=pnts.begin();
73			it!=pnts.end(); ++it) {
74	auterman	1.1	if ( func( *it, mass) )
75			result->SetPoint(i++, it->MZERO, it->MHALF);
76			}
77			return result;
78			}
79
80	auterman	1.1.1.2	template <class T> TH2F * TScan::Area(std::vector<T> pnts, double(*func)(T))
81	auterman	1.1	{
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	auterman	1.1.1.2	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	auterman	1.1	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	auterman	1.1.1.2	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	auterman	1.1
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	auterman	1.1.1.2	///Read in the generated x-sections
321			ReadGeneratedXsects( "xsec_scan.dat" );
322			cout << "Read in " << xsects.size() << " x-section points." <<endl;
323
324	auterman	1.1	///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	auterman	1.1.1.2	TH2F * ch_LSP = Area<SUSY_POINT>( points, charged_LSP );
332	auterman	1.1	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	auterman	1.1.1.2	c1->SetRightMargin ( 0.15 );
341	auterman	1.1	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	auterman	1.1.1.2	//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	auterman	1.1	plane->GetYaxis()->SetTitleOffset( 1.8 );
353	auterman	1.1.1.2	plane->GetZaxis()->SetTitleOffset( 1.1 );
354			plane->GetZaxis()->SetRangeUser(0.00001, 1000.);
355			gPad->SetLogz(1);
356			plane->Draw("COLZ");
357	auterman	1.1
358			///Now, add what ever should be displayed else
359			//---------------------------------------------------------------------------//
360	auterman	1.1.1.2	TLatex t;
361			t.SetTextColor(14);
362			t.SetTextAlign(12);
363			t.SetTextSize(0.02);
364	auterman	1.1
365			///gluino
366	auterman	1.1.1.2	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	auterman	1.1
377			///neutralino
378			//MyTGraph * m_neutralino1 = IsoMassLine( neutralino_1, 500. ); m_neutralino1->Draw();
379
380			///left-handed up-like squarks
381	auterman	1.1.1.2	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	auterman	1.1
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	auterman	1.1.1.2	t.SetTextSize(0.03);t.SetTextColor( 6 );
412			t.DrawLatex(200.,1800.,"charged LSP (#tilde{#tau})");
413	auterman	1.1
414	auterman	1.1.1.2	// 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	auterman	1.1	c1->SaveAs("dirt.eps");
471	auterman	1.1.1.2	gApplication->Run(); //interactive
472
473	auterman	1.1	//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	auterman	1.1.1.2
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	auterman	1.1
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			}