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


int TScan::DoStuff()
{
   ///Read in the generated masses
   ReadGeneratedMasses( "mass_scan.dat" );
   cout << "Read in " << points.size() << " mass 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();


   //c2->Draw();             //interactive
   //gApplication->Run();    //interactive

   c1->SaveAs("dirt.eps");
   //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.1.1.1 (vendor branch)
Committed:	Fri Nov 30 13:10:50 2007 UTC (17 years, 5 months ago) by auterman
Content type:	text/plain
Branch:	tex, Demo, SusyScan
Changes since 1.1:	+0 -0 lines
Log Message:	Scripts for a SUSY parameter scan
#	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
180	int TScan::DoStuff()
181	{
182	///Read in the generated masses
183	ReadGeneratedMasses( "mass_scan.dat" );
184	cout << "Read in " << points.size() << " mass points." <<endl;
185
186	///Where are no mass-points generated?
187	TH2F * no_solution = Area51( );
188	MyTGraph * contour_no_solution = GetContour( no_solution );
189	contour_no_solution->SetLineColor(1);
190	contour_no_solution->SetLineWidth(4);
191
192	///Where is the LSP charged?
193	TH2F * ch_LSP = Area( charged_LSP );
194	MyTGraph * contour_charged_LSP = GetContour( ch_LSP );
195	contour_charged_LSP->SetLineColor(6);
196	contour_charged_LSP->SetLineWidth(4);
197
198	///Main Canvas
199	TCanvas * c1 = new TCanvas("c1","c1",600,600);
200	c1->SetFillStyle ( 4000 );
201	c1->SetLeftMargin ( 0.15 );
202	c1->SetRightMargin ( 0.05 );
203	c1->SetBottomMargin( 0.10 );
204	c1->cd();
205	//c1->SetTopMargin ( 0.05 );
206	//TPostScript ps("plots.ps",-112);
207
208	///The base hist of the mSUGRA-plane with axis labels
209	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);
210	plane->GetYaxis()->SetTitleOffset( 1.8 );
211	plane->Draw();
212
213	///Now, add what ever should be displayed else
214	//---------------------------------------------------------------------------//
215
216	///gluino
217	//MyTGraph * m_gluino = IsoMassLine( gluino, 1000. ); m_gluino->Draw();
218
219	///neutralino
220	//MyTGraph * m_neutralino1 = IsoMassLine( neutralino_1, 500. ); m_neutralino1->Draw();
221
222	///left-handed up-like squarks
223	MyTGraph * m_ul5 = IsoMassLine( ul, 500. ); m_ul5->Draw();
224	MyTGraph * m_ul10 = IsoMassLine( ul, 1000. ); m_ul10->Draw();
225	MyTGraph * m_ul15 = IsoMassLine( ul, 1500. ); m_ul15->Draw();
226	MyTGraph * m_ul20 = IsoMassLine( ul, 2000. ); m_ul20->Draw();
227	MyTGraph * m_ul25 = IsoMassLine( ul, 2500. ); m_ul25->Draw();
228	MyTGraph * m_ul30 = IsoMassLine( ul, 3000. ); m_ul30->Draw();
229	MyTGraph * m_ul35 = IsoMassLine( ul, 3500. ); m_ul35->Draw();
230	MyTGraph * m_ul40 = IsoMassLine( ul, 4000. ); m_ul40->Draw();
231
232	///Contours
233	contour_no_solution->Add(max_x+5,430);
234	contour_no_solution->Add(max_x+5,min_y-5);
235	contour_no_solution->Add(min_x,min_y-5);
236	contour_no_solution->SetFillColor( 1 );
237	contour_no_solution->SetFillStyle( 3454 );
238	contour_no_solution->Draw("lf");
239	contour_no_solution->Draw();
240
241	contour_charged_LSP->Add(min_x-5, 726.25);
242	contour_charged_LSP->Add(min_x-5, max_y+5);
243	contour_charged_LSP->SetFillColor( 6 );
244	contour_charged_LSP->SetFillStyle( 3454 );
245	contour_charged_LSP->SetLineColor( 6 );
246	contour_charged_LSP->SetLineWidth( 4 );
247	contour_charged_LSP->Draw("lf");
248	contour_charged_LSP->Draw();
249
250
251
252	//c2->Draw(); //interactive
253	//gApplication->Run(); //interactive
254
255	c1->SaveAs("dirt.eps");
256	//ps.NewPage();
257
258	}
259
260	int scan(void)
261	{
262	gStyle->SetHistFillColor(0);
263	gStyle->SetPalette(1);
264	//gStyle->SetFillColor(0);
265	//gStyle->SetOptStat(kFALSE);
266	gStyle->SetCanvasColor(0);
267	gStyle->SetCanvasBorderMode(0);
268	gStyle->SetPadColor(0);
269	gStyle->SetPadBorderMode(0);
270	gStyle->SetFrameBorderMode(0);
271
272	gStyle->SetTitleFillColor(0);
273	gStyle->SetTitleBorderSize(0);
274	gStyle->SetTitleX(0.10);
275	gStyle->SetTitleY(0.98);
276	gStyle->SetTitleW(0.8);
277	gStyle->SetTitleH(0.06);
278
279	gStyle->SetErrorX(0);
280	gStyle->SetStatColor(0);
281	gStyle->SetStatBorderSize(0);
282	gStyle->SetStatX(0);
283	gStyle->SetStatY(0);
284	gStyle->SetStatW(0);
285	gStyle->SetStatH(0);
286
287	TScan * scan = new TScan();
288	scan->DoStuff();
289	delete scan;
290	}
291
292	int main(int argc, char** argv)
293	{
294	TApplication theApp("App", &argc, argv);
295	return scan();
296	}