SusyScan/PlotScript/plot.cc

#include "plot.h"
#include "SusyScan.h"
#include "GeneratorMasses.h"
#include "IsoMassLine.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->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");
   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 * doIt = new TheLimits();
   for (int i = 1; i<argc; ++i)
   {
     doIt->add( new SusyScan(argv[i]) );
   }
   doIt->match(genmasses);

   //c1->SetLogz(1);

   IsoMassLine<SusyScan> * mline = new IsoMassLine<SusyScan>(doIt->GetScan());

   TH2F*h = new TH2F("xsec",";m_{0} [GeV]; m_{1/2} [GeV]; x-section [pb]",
                     100,0,1009.9,50,0,500);
   //doIt->plot(h, Mzero, Mhalf, ObsXsecLimit);
   //doIt->plot(h, Mzero, Mhalf, Xsection);
   //doIt->plot(h, Mzero, Mhalf, XsecOverObserved);
   //doIt->plot(h, Mzero, Mhalf, XsecOverExpected);
   //doIt->plot(h, Mzero, Mhalf, SignalAcceptance);
   //doIt->plot(h, Mzero, Mhalf, NSignExpLimit);
   doIt->plot(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 = mline->get(Mzero, Mhalf, MGluino, 500);
   gl500->Draw();

   TGraph * sq500 = mline->get(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);
   doIt->plot(h_qg, MSquarkL, MGluino, SignalAcceptance);
   doIt->plot(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)
{
 // TRint *theApp = new TRint("ROOT example", &argc, argv);
  
  return plot(argc, argv);
}
Revision:	1.3
Committed:	Sat Nov 27 11:19:35 2010 UTC (14 years, 5 months ago) by auterman
Content type:	text/plain
Branch:	MAIN
Changes since 1.2:	+1 -1 lines
Log Message:	templating IsoMassLine; making specialized sort function private member
#	User	Rev	Content
1	auterman	1.1	#include "plot.h"
2			#include "SusyScan.h"
3	auterman	1.2	#include "GeneratorMasses.h"
4			#include "IsoMassLine.h"
5			#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->SetFillColor(0);
61			//gStyle->SetOptStat(kFALSE);
62			gStyle->SetCanvasColor(0);
63			gStyle->SetCanvasBorderMode(0);
64			gStyle->SetPadColor(0);
65			gStyle->SetPadBorderMode(0);
66			gStyle->SetFrameBorderMode(0);
67
68			gStyle->SetTitleFillColor(0);
69			gStyle->SetTitleBorderSize(0);
70			gStyle->SetTitleX(0.10);
71			gStyle->SetTitleY(0.98);
72			gStyle->SetTitleW(0.8);
73			gStyle->SetTitleH(0.06);
74
75			gStyle->SetErrorX(0);
76			gStyle->SetStatColor(0);
77			gStyle->SetStatBorderSize(0);
78			gStyle->SetStatX(0);
79			gStyle->SetStatY(0);
80			gStyle->SetStatW(0);
81			gStyle->SetStatH(0);
82
83			gStyle->SetTitleFont(22);
84			gStyle->SetLabelFont(22,"X");
85			gStyle->SetLabelFont(22,"Y");
86			gStyle->SetLabelFont(22,"Z");
87			gStyle->SetLabelSize(0.03,"X");
88			gStyle->SetLabelSize(0.03,"Y");
89			gStyle->SetLabelSize(0.03,"Z");
90
91			//gROOT->SetStyle("MyStyle");
92			TCanvas * c1 = new TCanvas("c1","c1",600,600);
93			c1->SetFillStyle ( 4000 );
94			c1->SetLeftMargin ( 0.15 );
95			c1->SetRightMargin ( 0.15 );
96			c1->SetBottomMargin( 0.10 );
97			c1->cd();
98
99	auterman	1.2	std::vector<GeneratorMasses> genmasses = FillGeneratorMasses("tb10_mu1_a0_massscan.dat");
100	auterman	1.1	TheLimits * doIt = new TheLimits();
101			for (int i = 1; i<argc; ++i)
102			{
103			doIt->add( new SusyScan(argv[i]) );
104			}
105	auterman	1.2	doIt->match(genmasses);
106
107			//c1->SetLogz(1);
108
109	auterman	1.3	IsoMassLine<SusyScan> * mline = new IsoMassLine<SusyScan>(doIt->GetScan());
110	auterman	1.2
111			TH2F*h = new TH2F("xsec",";m_{0} [GeV]; m_{1/2} [GeV]; x-section [pb]",
112			100,0,1009.9,50,0,500);
113			//doIt->plot(h, Mzero, Mhalf, ObsXsecLimit);
114			//doIt->plot(h, Mzero, Mhalf, Xsection);
115			//doIt->plot(h, Mzero, Mhalf, XsecOverObserved);
116			//doIt->plot(h, Mzero, Mhalf, XsecOverExpected);
117			//doIt->plot(h, Mzero, Mhalf, SignalAcceptance);
118			//doIt->plot(h, Mzero, Mhalf, NSignExpLimit);
119			doIt->plot(h, Mzero, Mhalf, ObsExclusion);
120			//h->SetMaximum(27);
121			//h->SetMinimum(0.01);
122			//c1->SetLogz(1);
123			h->GetYaxis()->SetTitleOffset(1.3);
124			h->GetZaxis()->SetTitleOffset(1.2);
125			h->Draw("colz");
126	auterman	1.1
127	auterman	1.2	TGraph * gl500 = mline->get(Mzero, Mhalf, MGluino, 500);
128			gl500->Draw();
129	auterman	1.1
130	auterman	1.2	TGraph * sq500 = mline->get(Mzero, Mhalf, MSquarkL, 500, 10);
131			sq500->Draw();
132
133			TH2F*h_qg = new TH2F("AccMGMSQ",";m_{#tilde{q}} [GeV]; m_{#tilde{g}} [GeV]; signal acceptance",
134			60,200,1400,50,200,1200);
135			doIt->plot(h_qg, MSquarkL, MGluino, SignalAcceptance);
136			doIt->plot(h_qg, MSquarkL, MGluino, NSignExpLimit);
137			//h_qg->SetMaximum(1.0);
138			h_qg->SetMinimum(0.01);
139
140			//c1->SetLogz(1);
141			h_qg->GetYaxis()->SetTitleOffset(1.3);
142			//h_qg->Draw("colz");
143
144
145	auterman	1.1
146	auterman	1.2	//c1->SaveAs("plot.pdf");
147
148
149			theApp.Run();
150	auterman	1.1	}
151
152
153			int main(int argc, char** argv)
154			{
155			// TRint *theApp = new TRint("ROOT example", &argc, argv);
156
157			return plot(argc, argv);
158			}