SusyScan/PlotScript/plot.cc

#include "plot.h"
#include "SusyScan.h"
#include "GeneratorMasses.h"
#include "PlotTools.h"
#include "TheLimits.h"

#include "TRint.h"
#include "TROOT.h"
#include "TObjArray.h"
#include "TStyle.h"

#include "TChain.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH1.h"
#include "TH2.h"
#include "TH2F.h"
#include "TTree.h"
#include "TKey.h"
#include "Riostream.h"
#include "TCanvas.h"

#include <string>
#include <cmath>

const double Luminosity = 31.2; //[pb^-1]
double Mzero(const SusyScan* p){ return p->Mzero; }
double Mhalf(const SusyScan* p){ return p->Mhalf; }
double MGluino(const SusyScan* p){ return p->MGL; }
double MSquarkL(const SusyScan* p){ return p->MUL; }
double MSquarkR(const SusyScan* p){ return p->MUR; }
double MChi1(const SusyScan* p){ return p->MZ1; }
double MChi2(const SusyScan* p){ return p->MZ2; }
double MChi3(const SusyScan* p){ return p->MZ3; }
double MChi4(const SusyScan* p){ return p->MZ4; }
double MCha1(const SusyScan* p){ return p->MW1; }
double MCha2(const SusyScan* p){ return p->MW2; }
double Xsection(const SusyScan* p){ return p->Xsection; }
double ExpXsecLimit(const SusyScan* p){ return p->ExpXsecLimit; }
double ObsXsecLimit(const SusyScan* p){ return p->ObsXsecLimit; }
double ExpExclusion(const SusyScan* p){ return (ExpXsecLimit(p)<Xsection(p)&&ExpXsecLimit(p)>0.01?1:-1); }
double ObsExclusion(const SusyScan* p){ return (ObsXsecLimit(p)<Xsection(p)&&ObsXsecLimit(p)>0.01?1:-1); }
double SoverSqrtB(const SusyScan* p){ return p->signal/(sqrt(p->background)+p->background_uncertainty+p->signal_uncertainty); }
double XsecOverObserved(const SusyScan* p){ return (ObsXsecLimit(p)==0 ? 9999. : Xsection(p)/ObsXsecLimit(p)); }
double XsecOverExpected(const SusyScan* p){ return (ObsXsecLimit(p)==0 ? 9999. : Xsection(p)/ObsXsecLimit(p)); }
double SignalAcceptance(const SusyScan* p){ return  p->signal / (Luminosity*Xsection(p)); }
double NSignExpLimit(const SusyScan* p){ return  p->signal * ExpXsecLimit(p)/Xsection(p); }
  

int plot(int argc, char** argv)
{
   TApplication theApp("App", 0, 0);
   if (gROOT->IsBatch()) {
      fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
      return 1;
   }   

   gStyle->SetHistFillColor(0);
   gStyle->SetPalette(1);
   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");
   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();
   
   std::vector<GeneratorMasses> genmasses = FillGeneratorMasses("tb10_mu1_a0_massscan.dat");
   TheLimits * genpoints = new TheLimits();
   for (int i = 1; i<argc; ++i)
   {
     genpoints->add( new SusyScan(argv[i]) );
   }
   genpoints->match(genmasses);

   //c1->SetLogz(1);

   PlotTools<SusyScan> * plotTools = new PlotTools<SusyScan>(genpoints->GetScan());

   TH2F*h = new TH2F("xsec",";m_{0} [GeV]; m_{1/2} [GeV]; x-section [pb]",
                     100,0,1009.9,50,0,500);
   //plotTools->Area(h, Mzero, Mhalf, ObsXsecLimit);
   //plotTools->Area(h, Mzero, Mhalf, Xsection);
   //plotTools->Area(h, Mzero, Mhalf, XsecOverObserved);
   //plotTools->Area(h, Mzero, Mhalf, XsecOverExpected);
   //plotTools->Area(h, Mzero, Mhalf, SignalAcceptance);
   //plotTools->Area(h, Mzero, Mhalf, NSignExpLimit);
   plotTools->Area(h, Mzero, Mhalf, ObsExclusion);
   //h->SetMaximum(27);
   //h->SetMinimum(0.01);
   //c1->SetLogz(1);
   h->GetYaxis()->SetTitleOffset(1.3);
   h->GetZaxis()->SetTitleOffset(1.2);
   h->Draw("colz");
   
   TGraph * gl500 = plotTools->Line(Mzero, Mhalf, MGluino, 500);
   gl500->Draw();

   TGraph * sq500 = plotTools->Line(Mzero, Mhalf, MSquarkL, 500, 10);
   sq500->Draw();

   TH2F*h_qg = new TH2F("AccMGMSQ",";m_{#tilde{q}} [GeV]; m_{#tilde{g}} [GeV]; signal acceptance",
                     60,200,1400,50,200,1200);
   plotTools->Area(h_qg, MSquarkL, MGluino, SignalAcceptance);
   plotTools->Area(h_qg, MSquarkL, MGluino, NSignExpLimit);
   //h_qg->SetMaximum(1.0);
   h_qg->SetMinimum(0.01);

   //c1->SetLogz(1);
   h_qg->GetYaxis()->SetTitleOffset(1.3);
   //h_qg->Draw("colz");
   
   
   //c1->SaveAs("plot.pdf");


   theApp.Run();
}


int main(int argc, char** argv)
{
  return plot(argc, argv);
}
Revision:	1.4
Committed:	Sat Nov 27 12:03:26 2010 UTC (14 years, 5 months ago) by auterman
Content type:	text/plain
Branch:	MAIN
Changes since 1.3:	+16 -20 lines
Log Message:	reorganization of code
#	User	Rev	Content
1	auterman	1.1	#include "plot.h"
2			#include "SusyScan.h"
3	auterman	1.2	#include "GeneratorMasses.h"
4	auterman	1.4	#include "PlotTools.h"
5	auterman	1.2	#include "TheLimits.h"
6	auterman	1.1
7			#include "TRint.h"
8			#include "TROOT.h"
9			#include "TObjArray.h"
10			#include "TStyle.h"
11
12			#include "TChain.h"
13			#include "TFile.h"
14	auterman	1.2	#include "TGraph.h"
15	auterman	1.1	#include "TH1.h"
16			#include "TH2.h"
17			#include "TH2F.h"
18			#include "TTree.h"
19			#include "TKey.h"
20			#include "Riostream.h"
21			#include "TCanvas.h"
22
23			#include <string>
24	auterman	1.2	#include <cmath>
25	auterman	1.1
26	auterman	1.2	const double Luminosity = 31.2; //[pb^-1]
27	auterman	1.1	double Mzero(const SusyScan* p){ return p->Mzero; }
28			double Mhalf(const SusyScan* p){ return p->Mhalf; }
29	auterman	1.2	double MGluino(const SusyScan* p){ return p->MGL; }
30			double MSquarkL(const SusyScan* p){ return p->MUL; }
31			double MSquarkR(const SusyScan* p){ return p->MUR; }
32			double MChi1(const SusyScan* p){ return p->MZ1; }
33			double MChi2(const SusyScan* p){ return p->MZ2; }
34			double MChi3(const SusyScan* p){ return p->MZ3; }
35			double MChi4(const SusyScan* p){ return p->MZ4; }
36			double MCha1(const SusyScan* p){ return p->MW1; }
37			double MCha2(const SusyScan* p){ return p->MW2; }
38	auterman	1.1	double Xsection(const SusyScan* p){ return p->Xsection; }
39	auterman	1.2	double ExpXsecLimit(const SusyScan* p){ return p->ExpXsecLimit; }
40			double ObsXsecLimit(const SusyScan* p){ return p->ObsXsecLimit; }
41			double ExpExclusion(const SusyScan* p){ return (ExpXsecLimit(p)<Xsection(p)&&ExpXsecLimit(p)>0.01?1:-1); }
42			double ObsExclusion(const SusyScan* p){ return (ObsXsecLimit(p)<Xsection(p)&&ObsXsecLimit(p)>0.01?1:-1); }
43			double SoverSqrtB(const SusyScan* p){ return p->signal/(sqrt(p->background)+p->background_uncertainty+p->signal_uncertainty); }
44			double XsecOverObserved(const SusyScan* p){ return (ObsXsecLimit(p)==0 ? 9999. : Xsection(p)/ObsXsecLimit(p)); }
45			double XsecOverExpected(const SusyScan* p){ return (ObsXsecLimit(p)==0 ? 9999. : Xsection(p)/ObsXsecLimit(p)); }
46			double SignalAcceptance(const SusyScan* p){ return p->signal / (Luminosity*Xsection(p)); }
47			double NSignExpLimit(const SusyScan* p){ return p->signal * ExpXsecLimit(p)/Xsection(p); }
48	auterman	1.1
49
50			int plot(int argc, char** argv)
51			{
52	auterman	1.2	TApplication theApp("App", 0, 0);
53			if (gROOT->IsBatch()) {
54			fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
55			return 1;
56			}
57
58	auterman	1.1	gStyle->SetHistFillColor(0);
59			gStyle->SetPalette(1);
60			gStyle->SetCanvasColor(0);
61			gStyle->SetCanvasBorderMode(0);
62			gStyle->SetPadColor(0);
63			gStyle->SetPadBorderMode(0);
64			gStyle->SetFrameBorderMode(0);
65
66			gStyle->SetTitleFillColor(0);
67			gStyle->SetTitleBorderSize(0);
68			gStyle->SetTitleX(0.10);
69			gStyle->SetTitleY(0.98);
70			gStyle->SetTitleW(0.8);
71			gStyle->SetTitleH(0.06);
72
73			gStyle->SetErrorX(0);
74			gStyle->SetStatColor(0);
75			gStyle->SetStatBorderSize(0);
76			gStyle->SetStatX(0);
77			gStyle->SetStatY(0);
78			gStyle->SetStatW(0);
79			gStyle->SetStatH(0);
80
81			gStyle->SetTitleFont(22);
82			gStyle->SetLabelFont(22,"X");
83			gStyle->SetLabelFont(22,"Y");
84			gStyle->SetLabelFont(22,"Z");
85			gStyle->SetLabelSize(0.03,"X");
86			gStyle->SetLabelSize(0.03,"Y");
87			gStyle->SetLabelSize(0.03,"Z");
88
89			//gROOT->SetStyle("MyStyle");
90			TCanvas * c1 = new TCanvas("c1","c1",600,600);
91			c1->SetFillStyle ( 4000 );
92			c1->SetLeftMargin ( 0.15 );
93			c1->SetRightMargin ( 0.15 );
94			c1->SetBottomMargin( 0.10 );
95			c1->cd();
96
97	auterman	1.2	std::vector<GeneratorMasses> genmasses = FillGeneratorMasses("tb10_mu1_a0_massscan.dat");
98	auterman	1.4	TheLimits * genpoints = new TheLimits();
99	auterman	1.1	for (int i = 1; i<argc; ++i)
100			{
101	auterman	1.4	genpoints->add( new SusyScan(argv[i]) );
102	auterman	1.1	}
103	auterman	1.4	genpoints->match(genmasses);
104	auterman	1.2
105			//c1->SetLogz(1);
106
107	auterman	1.4	PlotTools<SusyScan> * plotTools = new PlotTools<SusyScan>(genpoints->GetScan());
108	auterman	1.2
109			TH2F*h = new TH2F("xsec",";m_{0} [GeV]; m_{1/2} [GeV]; x-section [pb]",
110			100,0,1009.9,50,0,500);
111	auterman	1.4	//plotTools->Area(h, Mzero, Mhalf, ObsXsecLimit);
112			//plotTools->Area(h, Mzero, Mhalf, Xsection);
113			//plotTools->Area(h, Mzero, Mhalf, XsecOverObserved);
114			//plotTools->Area(h, Mzero, Mhalf, XsecOverExpected);
115			//plotTools->Area(h, Mzero, Mhalf, SignalAcceptance);
116			//plotTools->Area(h, Mzero, Mhalf, NSignExpLimit);
117			plotTools->Area(h, Mzero, Mhalf, ObsExclusion);
118	auterman	1.2	//h->SetMaximum(27);
119			//h->SetMinimum(0.01);
120			//c1->SetLogz(1);
121			h->GetYaxis()->SetTitleOffset(1.3);
122			h->GetZaxis()->SetTitleOffset(1.2);
123			h->Draw("colz");
124	auterman	1.1
125	auterman	1.4	TGraph * gl500 = plotTools->Line(Mzero, Mhalf, MGluino, 500);
126	auterman	1.2	gl500->Draw();
127	auterman	1.1
128	auterman	1.4	TGraph * sq500 = plotTools->Line(Mzero, Mhalf, MSquarkL, 500, 10);
129	auterman	1.2	sq500->Draw();
130
131			TH2F*h_qg = new TH2F("AccMGMSQ",";m_{#tilde{q}} [GeV]; m_{#tilde{g}} [GeV]; signal acceptance",
132			60,200,1400,50,200,1200);
133	auterman	1.4	plotTools->Area(h_qg, MSquarkL, MGluino, SignalAcceptance);
134			plotTools->Area(h_qg, MSquarkL, MGluino, NSignExpLimit);
135	auterman	1.2	//h_qg->SetMaximum(1.0);
136			h_qg->SetMinimum(0.01);
137
138			//c1->SetLogz(1);
139			h_qg->GetYaxis()->SetTitleOffset(1.3);
140			//h_qg->Draw("colz");
141
142
143	auterman	1.1
144	auterman	1.2	//c1->SaveAs("plot.pdf");
145
146
147			theApp.Run();
148	auterman	1.1	}
149
150
151			int main(int argc, char** argv)
152			{
153			return plot(argc, argv);
154			}