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 "scan.h"


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

using namespace std;

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)<1.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 ); }


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()
{
}

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

TH2F * TScan::Area(bool(*func)(SUSY_POINT))
{
  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 (std::vector<SUSY_POINT>::const_iterator it=points.begin();
       it!=points.end(); ++it) {
     if ( func( *it ) ) 
       result->Fill(it->MZERO, it->MHALF, 1.);
//if (it->MZERO==100. && it->MHALF==200.)  
//cout << it->MZERO <<", "<<it->MHALF
//     << ",  m(Z1)="<<it->MZ1<<", m(stau)="<<it->MTAU1<<endl;
  } 
  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( 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.05 );
   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);
   plane->GetYaxis()->SetTitleOffset(  1.8 );
   plane->Draw();

   ///Now, add what ever should be displayed else
   //---------------------------------------------------------------------------//
   
   ///gluino
   //MyTGraph * m_gluino = IsoMassLine( gluino, 1000. ); m_gluino->Draw();
   
   ///neutralino
   //MyTGraph * m_neutralino1 = IsoMassLine( neutralino_1, 500. ); m_neutralino1->Draw();
   
   ///left-handed up-like squarks
   MyTGraph * m_ul5  = IsoMassLine( ul,  500. ); m_ul5->Draw();
   MyTGraph * m_ul10 = IsoMassLine( ul, 1000. ); m_ul10->Draw();
   MyTGraph * m_ul15 = IsoMassLine( ul, 1500. ); m_ul15->Draw();
   MyTGraph * m_ul20 = IsoMassLine( ul, 2000. ); m_ul20->Draw();
   MyTGraph * m_ul25 = IsoMassLine( ul, 2500. ); m_ul25->Draw();
   MyTGraph * m_ul30 = IsoMassLine( ul, 3000. ); m_ul30->Draw();
   MyTGraph * m_ul35 = IsoMassLine( ul, 3500. ); m_ul35->Draw();
   MyTGraph * m_ul40 = IsoMassLine( ul, 4000. ); m_ul40->Draw();

   ///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();


   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);
   
   TScan * scan = new TScan();
   scan->DoStuff();
   delete scan;
}

int main(int argc, char** argv)
{
  TApplication theApp("App", &argc, argv);
  return scan();
}
Revision:	1.2
Committed:	Tue Dec 4 12:57:34 2007 UTC (17 years, 5 months ago) by auterman
Content type:	text/plain
Branch:	MAIN
Changes since 1.1:	+116 -3 lines
Log Message:	* empty log message *
#	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
18
19	#include "scan.h"
20
21
22	TROOT root("GUI", "GUI test environement");
23
24	using namespace std;
25
26	bool charged_LSP(const SUSY_POINT point)
27	{ return ( (fabs(point.MZ1) - point.MTAU1)>0.0 ); }
28
29	bool gluino(const SUSY_POINT point, const double mass)
30	{ return ( fabs(point.MGL - mass)<1.0 ); }
31
32	bool ul(const SUSY_POINT point, const double mass)
33	{ return ( fabs(point.MUL - mass)<5.0 ); }
34
35	bool neutralino_1(const SUSY_POINT point, const double mass)
36	{ return ( fabs(fabs(point.MZ1) - mass)<1.0 ); }
37
38
39	TScan::TScan()
40	{
41	plot_id = 0;
42	bins_x = 200;
43	min_x = 0.0;
44	max_x = 2000.0;
45	bins_y = 200;
46	min_x = 0.0;
47	max_y = 2000.0;
48	}
49
50	TScan::~TScan()
51	{
52	}
53
54	MyTGraph * TScan::IsoMassLine(bool(*func)(SUSY_POINT,double), double mass )
55	{
56	MyTGraph * result = new MyTGraph(1);
57	int i=0;
58	for (std::vector<SUSY_POINT>::const_iterator it=points.begin();
59	it!=points.end(); ++it) {
60	if ( func( *it, mass) )
61	result->SetPoint(i++, it->MZERO, it->MHALF);
62	}
63	return result;
64	}
65
66	TH2F * TScan::Area(bool(*func)(SUSY_POINT))
67	{
68	char * name = new char[255];
69	sprintf(name,"Area_plot%d",++plot_id);
70	TH2F * result = new TH2F(name,"",bins_x,min_x,max_x,bins_y,min_x,max_y);
71	for (std::vector<SUSY_POINT>::const_iterator it=points.begin();
72	it!=points.end(); ++it) {
73	if ( func( *it ) )
74	result->Fill(it->MZERO, it->MHALF, 1.);
75	//if (it->MZERO==100. && it->MHALF==200.)
76	//cout << it->MZERO <<", "<<it->MHALF
77	// << ", m(Z1)="<<it->MZ1<<", m(stau)="<<it->MTAU1<<endl;
78	}
79	delete name;
80	return result;
81	}
82
83	TH2F * TScan::Area51()
84	{
85	char * name = new char[255];
86	sprintf(name,"Area51_plot%d",++plot_id);
87	TH2F * result = new TH2F(name,"",bins_x,min_x,max_x,bins_y,min_x,max_y);
88	for (int x=1; x<=bins_x; ++x){
89	for (int y=1; y<=bins_y; ++y){
90	result->SetBinContent(x,y,1.);
91	}
92	}
93	for (std::vector<SUSY_POINT>::const_iterator it=points.begin();
94	it!=points.end(); ++it) {
95	result->SetBinContent((int)(bins_x*it->MZERO/(max_x-min_x)),
96	(int)(bins_y*it->MHALF/(max_y-min_y)),0.);
97	}
98	delete name;
99	return result;
100	}
101
102	MyTGraph * TScan::GetContour(TH2F*plot, int flag)
103	{
104	char * name = new char[255];
105	sprintf(name,"canv%d",++plot_id);
106	TCanvas * canv = new TCanvas(name,"c1",600,600);
107	canv->cd();
108	plot->SetContour(2);
109	plot->SetContourLevel(0,0);
110	plot->SetContourLevel(1,0.5);
111	plot->SetFillColor(1);
112	plot->Draw("CONT List");
113	canv->Update();
114	TObjArray contours = (TObjArray)gROOT->GetListOfSpecials()->FindObject("contours");
115	int ncontours = contours->GetSize();
116	if (ncontours==0) return 0;
117	//cout << "N contours = " << ncontours << endl;
118	TList list = (TList)contours->At(0);
119	MyTGraph gr1 = (MyTGraph)list->First();
120	delete name;
121	delete canv;
122	delete list;
123
124	return gr1;
125	}
126
127	void TScan::ReadGeneratedMasses(std::string file)
128	{
129	ifstream masses_file;
130	masses_file.open(file.c_str());
131	SUSY_POINT p;
132	while (1) {
133	masses_file >> p.MZERO
134	>> p.MHALF
135	>> p.TANB
136	>> p.SGNMU
137	>> p.AZERO
138	>> p.MTOP
139	>> p.muQ
140	>> p.Q
141	>> p.M1
142	>> p.M2
143	>> p.M3
144	>> p.MGL
145	>> p.MUL
146	>> p.MB1
147	>> p.MSN
148	>> p.MNTAU
149	>> p.MZ1
150	>> p.MW1
151	>> p.MHL
152	>> p.MUR
153	>> p.MB2
154	>> p.MEL
155	>> p.MTAU1
156	>> p.MZ2
157	>> p.MW2
158	>> p.MHH
159	>> p.MDL
160	>> p.MT1
161	>> p.MER
162	>> p.MTAU2
163	>> p.MZ3
164	>> p.MHA
165	>> p.MDR
166	>> p.MT2
167	>> p.MZ4
168	>> p.MHp;
169
170	if (!masses_file.good()) break;
171	if (fabs(p.SGNMU)!=1.) {
172	cerr << "check lines near m0=" << p.MZERO << ", m1/2=" << p.MHALF << endl;
173	break;
174	}
175	points.push_back(p);
176	}
177	}
178
179	void TScan::ReadGeneratedXsects(std::string file)
180	{
181	ifstream xsects_file;
182	xsects_file.open(file.c_str());
183	SUSY_XSECS x;
184	while (1) {
185	xsects_file >> x.MZERO
186	>> x.MHALF
187	>> x.TANB
188	>> x.SGNMU
189	>> x.AZERO
190	>> x.MTOP
191	>> x.XS0
192	>> x.XS201
193	>> x.XS202
194	>> x.XS204
195	>> x.XS205
196	>> x.XS207
197	>> x.XS208
198	>> x.XS209
199	>> x.XS210
200	>> x.XS211
201	>> x.XS212
202	>> x.XS213
203	>> x.XS214
204	>> x.XS216
205	>> x.XS217
206	>> x.XS218
207	>> x.XS219
208	>> x.XS220
209	>> x.XS221
210	>> x.XS222
211	>> x.XS223
212	>> x.XS224
213	>> x.XS225
214	>> x.XS226
215	>> x.XS227
216	>> x.XS228
217	>> x.XS229
218	>> x.XS230
219	>> x.XS231
220	>> x.XS232
221	>> x.XS233
222	>> x.XS234
223	>> x.XS235
224	>> x.XS236
225	>> x.XS237
226	>> x.XS238
227	>> x.XS239
228	>> x.XS240
229	>> x.XS241
230	>> x.XS242
231	>> x.XS243
232	>> x.XS244
233	>> x.XS246
234	>> x.XS247
235	>> x.XS248
236	>> x.XS249
237	>> x.XS250
238	>> x.XS251
239	>> x.XS252
240	>> x.XS253
241	>> x.XS254
242	>> x.XS256
243	>> x.XS258
244	>> x.XS259
245	>> x.XS261
246	>> x.XS262
247	>> x.XS263
248	>> x.XS264
249	>> x.XS265
250	>> x.XS271
251	>> x.XS272
252	>> x.XS273
253	>> x.XS274
254	>> x.XS275
255	>> x.XS276
256	>> x.XS277
257	>> x.XS278
258	>> x.XS279
259	>> x.XS280
260	>> x.XS281
261	>> x.XS282
262	>> x.XS283
263	>> x.XS284
264	>> x.XS285
265	>> x.XS286
266	>> x.XS287
267	>> x.XS288
268	>> x.XS289
269	>> x.XS290
270	>> x.XS291
271	>> x.XS292
272	>> x.XS293
273	>> x.XS294
274	>> x.XS295
275	>> x.XS296
276	>> x.XS297
277	>> x.XS298
278	>> x.XS299
279	>> x.XS300
280	>> x.XS301;
281	if (!xsects_file.good()) break;
282	if (fabs(x.SGNMU)!=1.) {
283	cerr << "check lines near m0=" << x.MZERO << ", m1/2=" << x.MHALF << endl;
284	break;
285	}
286	xsects.push_back(x);
287	}
288	}
289
290	int TScan::DoStuff()
291	{
292	///Read in the generated masses
293	ReadGeneratedMasses( "mass_scan.dat" );
294	cout << "Read in " << points.size() << " mass points." <<endl;
295
296	///Read in the generated x-sections
297	ReadGeneratedXsects( "xsec_scan.dat" );
298	cout << "Read in " << xsects.size() << " x-section points." <<endl;
299
300	///Where are no mass-points generated?
301	TH2F * no_solution = Area51( );
302	MyTGraph * contour_no_solution = GetContour( no_solution );
303	contour_no_solution->SetLineColor(1);
304	contour_no_solution->SetLineWidth(4);
305
306	///Where is the LSP charged?
307	TH2F * ch_LSP = Area( charged_LSP );
308	MyTGraph * contour_charged_LSP = GetContour( ch_LSP );
309	contour_charged_LSP->SetLineColor(6);
310	contour_charged_LSP->SetLineWidth(4);
311
312	///Main Canvas
313	TCanvas * c1 = new TCanvas("c1","c1",600,600);
314	c1->SetFillStyle ( 4000 );
315	c1->SetLeftMargin ( 0.15 );
316	c1->SetRightMargin ( 0.05 );
317	c1->SetBottomMargin( 0.10 );
318	c1->cd();
319	//c1->SetTopMargin ( 0.05 );
320	//TPostScript ps("plots.ps",-112);
321
322	///The base hist of the mSUGRA-plane with axis labels
323	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);
324	plane->GetYaxis()->SetTitleOffset( 1.8 );
325	plane->Draw();
326
327	///Now, add what ever should be displayed else
328	//---------------------------------------------------------------------------//
329
330	///gluino
331	//MyTGraph * m_gluino = IsoMassLine( gluino, 1000. ); m_gluino->Draw();
332
333	///neutralino
334	//MyTGraph * m_neutralino1 = IsoMassLine( neutralino_1, 500. ); m_neutralino1->Draw();
335
336	///left-handed up-like squarks
337	MyTGraph * m_ul5 = IsoMassLine( ul, 500. ); m_ul5->Draw();
338	MyTGraph * m_ul10 = IsoMassLine( ul, 1000. ); m_ul10->Draw();
339	MyTGraph * m_ul15 = IsoMassLine( ul, 1500. ); m_ul15->Draw();
340	MyTGraph * m_ul20 = IsoMassLine( ul, 2000. ); m_ul20->Draw();
341	MyTGraph * m_ul25 = IsoMassLine( ul, 2500. ); m_ul25->Draw();
342	MyTGraph * m_ul30 = IsoMassLine( ul, 3000. ); m_ul30->Draw();
343	MyTGraph * m_ul35 = IsoMassLine( ul, 3500. ); m_ul35->Draw();
344	MyTGraph * m_ul40 = IsoMassLine( ul, 4000. ); m_ul40->Draw();
345
346	///Contours
347	contour_no_solution->Add(max_x+5,430);
348	contour_no_solution->Add(max_x+5,min_y-5);
349	contour_no_solution->Add(min_x,min_y-5);
350	contour_no_solution->SetFillColor( 1 );
351	contour_no_solution->SetFillStyle( 3454 );
352	contour_no_solution->Draw("lf");
353	contour_no_solution->Draw();
354
355	contour_charged_LSP->Add(min_x-5, 726.25);
356	contour_charged_LSP->Add(min_x-5, max_y+5);
357	contour_charged_LSP->SetFillColor( 6 );
358	contour_charged_LSP->SetFillStyle( 3454 );
359	contour_charged_LSP->SetLineColor( 6 );
360	contour_charged_LSP->SetLineWidth( 4 );
361	contour_charged_LSP->Draw("lf");
362	contour_charged_LSP->Draw();
363
364
365
366	c1->SaveAs("dirt.eps");
367	gApplication->Run(); //interactive
368
369	//ps.NewPage();
370
371	}
372
373	int scan(void)
374	{
375	gStyle->SetHistFillColor(0);
376	gStyle->SetPalette(1);
377	//gStyle->SetFillColor(0);
378	//gStyle->SetOptStat(kFALSE);
379	gStyle->SetCanvasColor(0);
380	gStyle->SetCanvasBorderMode(0);
381	gStyle->SetPadColor(0);
382	gStyle->SetPadBorderMode(0);
383	gStyle->SetFrameBorderMode(0);
384
385	gStyle->SetTitleFillColor(0);
386	gStyle->SetTitleBorderSize(0);
387	gStyle->SetTitleX(0.10);
388	gStyle->SetTitleY(0.98);
389	gStyle->SetTitleW(0.8);
390	gStyle->SetTitleH(0.06);
391
392	gStyle->SetErrorX(0);
393	gStyle->SetStatColor(0);
394	gStyle->SetStatBorderSize(0);
395	gStyle->SetStatX(0);
396	gStyle->SetStatY(0);
397	gStyle->SetStatW(0);
398	gStyle->SetStatH(0);
399
400	TScan * scan = new TScan();
401	scan->DoStuff();
402	delete scan;
403	}
404
405	int main(int argc, char** argv)
406	{
407	TApplication theApp("App", &argc, argv);
408	return scan();
409	}