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
Committed:	Fri Nov 30 13:10:50 2007 UTC (17 years, 5 months ago) by auterman
Content type:	text/plain
Branch:	MAIN
Branch point for:	tex, Demo, SusyScan
Log Message:	Initial revision
#	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
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			}