ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/auterman/SusyScan/PlotScript/scan.cc
Revision: 1.1
Committed: Tue Nov 23 20:16:52 2010 UTC (14 years, 5 months ago) by auterman
Content type: text/plain
Branch: MAIN
Log Message:
adding new functionality

File Contents

# 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     #include <TLatex.h>
18     #include <TMarker.h>
19     #include "TRint.h"
20     #include "TROOT.h"
21     #include "TObjArray.h"
22    
23     #include "scan.h"
24     #include "SusyScan.h"
25    
26    
27     //TROOT root("GUI", "GUI test environement");
28    
29     using namespace std;
30    
31    
32     double total_xsec(const SUSY_XSECS point)
33     { return point.total(); }
34    
35    
36     bool charged_LSP(const SUSY_POINT point)
37     { return ( (fabs(point.MZ1) - point.MTAU1)>0.0 ); }
38    
39     bool gluino(const SUSY_POINT point, const double mass)
40     { return ( fabs(point.MGL - mass)<5.0 ); }
41    
42     bool ul(const SUSY_POINT point, const double mass)
43     { return ( fabs(point.MUL - mass)<5.0 ); }
44    
45     bool neutralino_1(const SUSY_POINT point, const double mass)
46     { return ( fabs(fabs(point.MZ1) - mass)<1.0 ); }
47    
48     bool sq_equal_gluino(const SUSY_POINT point, const double treshold)
49     { return ( fabs(point.MUL - point.MGL)<treshold ); }
50    
51    
52     bool xs_sq_equal_gluino(const SUSY_XSECS point, const double treshold)
53     { if (point.sgsg()<1E-5&&point.sqsq()<1E-5) return false; else return (fabs(point.sgsg() - point.sqsq())<treshold*point.sqsq()/100. ); }
54    
55    
56     TScan::TScan()
57     {
58     plot_id = 0;
59     bins_x = 200;
60     min_x = 0.0;
61     max_x = 2000.0;
62     bins_y = 40;
63     min_y = 0.0;
64     max_y = 1000.0;
65     }
66    
67     TScan::~TScan()
68     {
69     }
70    
71     template <class T> MyTGraph * TScan::IsoMassLine(std::vector<T> pnts, bool(*func)(T,double), double mass )
72     {
73     MyTGraph * result = new MyTGraph(1);
74     int i=0;
75     for (typename std::vector<T>::const_iterator it=pnts.begin();
76     it!=pnts.end(); ++it) {
77     if ( func( *it, mass) )
78     result->SetPoint(i++, it->MZERO, it->MHALF);
79     }
80     return result;
81     }
82    
83     template <class T> TH2F * TScan::Area(std::vector<T> pnts, double(*func)(T))
84     {
85     char * name = new char[255];
86     sprintf(name,"Area_plot%d",++plot_id);
87     TH2F * result = new TH2F(name,"",bins_x,min_x,max_x,bins_y,min_y,max_y);
88     for (typename std::vector<T>::const_iterator it=pnts.begin();
89     it!=pnts.end(); ++it) {
90     result->Fill(it->MZERO, it->MHALF, func(*it) );
91     }
92     delete name;
93     return result;
94     }
95    
96     template <class T> TH2F * TScan::Area(std::vector<T> pnts, bool(*func)(T))
97     {
98     char * name = new char[255];
99     sprintf(name,"Area_plot%d",++plot_id);
100     TH2F * result = new TH2F(name,"",bins_x,min_x,max_x,bins_y,min_y,max_y);
101     for (typename std::vector<T>::const_iterator it=pnts.begin();
102     it!=pnts.end(); ++it) {
103     if ( func( *it ) )
104     result->Fill(it->MZERO, it->MHALF, 1.);
105     }
106     delete name;
107     return result;
108     }
109    
110     TH2F * TScan::Area51()
111     {
112     char * name = new char[255];
113     sprintf(name,"Area51_plot%d",++plot_id);
114     TH2F * result = new TH2F(name,"",bins_x,min_x,max_x,bins_y,min_y,max_y);
115     for (int x=1; x<=bins_x; ++x){
116     for (int y=1; y<=bins_y; ++y){
117     result->SetBinContent(x,y,1.);
118     }
119     }
120     for (std::vector<SUSY_POINT>::const_iterator it=points.begin();
121     it!=points.end(); ++it) {
122     result->SetBinContent((int)(bins_x*it->MZERO/(max_x-min_x)),
123     (int)(bins_y*it->MHALF/(max_y-min_y)),0.);
124     }
125     delete name;
126     return result;
127     }
128    
129     TGraph * TScan::GetContour(TH2F*plot, int flag)
130     {
131     /*
132     char * name = new char[255];
133     sprintf(name,"canv%d",++plot_id);
134     TCanvas * canv = new TCanvas(name,"c1",600,600);
135     canv->cd();
136     plot->SetContour(2);
137     plot->SetContourLevel(0,0);
138     plot->SetContourLevel(1,0.5);
139     plot->SetFillColor(1);
140     plot->Draw("CONT List");
141     canv->Update();
142     TObjArray *contours = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
143     int ncontours = contours->GetSize();
144     if (ncontours==0) return 0;
145     //cout << "N contours = " << ncontours << endl;
146     TList *list = (TList*)contours->At(0);
147     MyTGraph *gr1 = (MyTGraph*)list->First();
148     delete name;
149     delete canv;
150     delete list;
151     */
152     plot->SetContour(2);
153     plot->SetContourLevel(0,0);
154     plot->SetContourLevel(1,0.5);
155     plot->SetFillColor(1);
156     plot->Draw("CONT List");
157     gPad->Update();
158     TObjArray *contours = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
159     int ncontours = contours->GetSize();
160     cout << "N contours = " << ncontours << endl;
161     TList *list = (TList*)contours->At(0);
162     //list->GetSize();
163     TGraph *gr1 = (TGraph*)list->First();
164     return gr1;
165     }
166    
167     void TScan::ReadGeneratedMasses(std::string file)
168     {
169     ifstream masses_file;
170     masses_file.open(file.c_str());
171     SUSY_POINT p;
172     while (1) {
173     masses_file >> p.MZERO
174     >> p.MHALF
175     >> p.TANB
176     >> p.SGNMU
177     >> p.AZERO
178     >> p.MTOP
179     >> p.muQ
180     >> p.Q
181     >> p.M1
182     >> p.M2
183     >> p.M3
184     >> p.MGL
185     >> p.MUL
186     >> p.MB1
187     >> p.MSN
188     >> p.MNTAU
189     >> p.MZ1
190     >> p.MW1
191     >> p.MHL
192     >> p.MUR
193     >> p.MB2
194     >> p.MEL
195     >> p.MTAU1
196     >> p.MZ2
197     >> p.MW2
198     >> p.MHH
199     >> p.MDL
200     >> p.MT1
201     >> p.MER
202     >> p.MTAU2
203     >> p.MZ3
204     >> p.MHA
205     >> p.MDR
206     >> p.MT2
207     >> p.MZ4
208     >> p.MHp;
209    
210     if (!masses_file.good()) break;
211     if (fabs(p.SGNMU)!=1.) {
212     cerr << "check lines near m0=" << p.MZERO << ", m1/2=" << p.MHALF << endl;
213     break;
214     }
215     points.push_back(p);
216     }
217     }
218    
219     void TScan::ReadGeneratedXsects(std::string file)
220     {
221     ifstream xsects_file;
222     xsects_file.open(file.c_str());
223     SUSY_XSECS x;
224     while (1) {
225     xsects_file >> x.MZERO
226     >> x.MHALF
227     >> x.TANB
228     >> x.SGNMU
229     >> x.AZERO
230     >> x.MTOP
231     >> x.XS0
232     >> x.XS201
233     >> x.XS202
234     >> x.XS204
235     >> x.XS205
236     >> x.XS207
237     >> x.XS208
238     >> x.XS209
239     >> x.XS210
240     >> x.XS211
241     >> x.XS212
242     >> x.XS213
243     >> x.XS214
244     >> x.XS216
245     >> x.XS217
246     >> x.XS218
247     >> x.XS219
248     >> x.XS220
249     >> x.XS221
250     >> x.XS222
251     >> x.XS223
252     >> x.XS224
253     >> x.XS225
254     >> x.XS226
255     >> x.XS227
256     >> x.XS228
257     >> x.XS229
258     >> x.XS230
259     >> x.XS231
260     >> x.XS232
261     >> x.XS233
262     >> x.XS234
263     >> x.XS235
264     >> x.XS236
265     >> x.XS237
266     >> x.XS238
267     >> x.XS239
268     >> x.XS240
269     >> x.XS241
270     >> x.XS242
271     >> x.XS243
272     >> x.XS244
273     >> x.XS246
274     >> x.XS247
275     >> x.XS248
276     >> x.XS249
277     >> x.XS250
278     >> x.XS251
279     >> x.XS252
280     >> x.XS253
281     >> x.XS254
282     >> x.XS256
283     >> x.XS258
284     >> x.XS259
285     >> x.XS261
286     >> x.XS262
287     >> x.XS263
288     >> x.XS264
289     >> x.XS265
290     >> x.XS271
291     >> x.XS272
292     >> x.XS273
293     >> x.XS274
294     >> x.XS275
295     >> x.XS276
296     >> x.XS277
297     >> x.XS278
298     >> x.XS279
299     >> x.XS280
300     >> x.XS281
301     >> x.XS282
302     >> x.XS283
303     >> x.XS284
304     >> x.XS285
305     >> x.XS286
306     >> x.XS287
307     >> x.XS288
308     >> x.XS289
309     >> x.XS290
310     >> x.XS291
311     >> x.XS292
312     >> x.XS293
313     >> x.XS294
314     >> x.XS295
315     >> x.XS296
316     >> x.XS297
317     >> x.XS298
318     >> x.XS299
319     >> x.XS300
320     >> x.XS301;
321     if (!xsects_file.good()) break;
322     if (fabs(x.SGNMU)!=1.) {
323     cerr << "check lines near m0=" << x.MZERO << ", m1/2=" << x.MHALF << endl;
324     break;
325     }
326     xsects.push_back(x);
327     }
328     }
329    
330     int TScan::DoStuff()
331     {
332     ///Read in the generated masses
333     ReadGeneratedMasses( "tb10_mu1_a0_massscan.dat" );
334     cout << "Read in " << points.size() << " mass points." <<endl;
335    
336     ///Read in the generated x-sections
337     ReadGeneratedXsects( "tb10mu1a0s7tev_scan_xsec.dat" );
338     cout << "Read in " << xsects.size() << " x-section points." <<endl;
339    
340     ///Where are no mass-points generated?
341     TH2F * no_solution = Area51( );
342     TGraph * contour_no_solution = GetContour( no_solution );
343     if (contour_no_solution) contour_no_solution->SetLineColor(1);
344     if (contour_no_solution) contour_no_solution->SetLineWidth(4);
345    
346     ///Where is the LSP charged?
347     TH2F * ch_LSP = Area<SUSY_POINT>( points, charged_LSP );
348     TGraph * contour_charged_LSP = GetContour( ch_LSP );
349     contour_charged_LSP->SetLineColor(6);
350     contour_charged_LSP->SetLineWidth(4);
351    
352     ///Main Canvas
353     TCanvas * c1 = new TCanvas("c1","c1",600,600);
354     c1->SetFillStyle ( 4000 );
355     c1->SetLeftMargin ( 0.15 );
356     c1->SetRightMargin ( 0.15 );
357     c1->SetBottomMargin( 0.10 );
358     c1->cd();
359     //c1->SetTopMargin ( 0.05 );
360     //TPostScript ps("plots.ps",-112);
361    
362     ///The base hist of the mSUGRA-plane with axis labels
363     //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);
364     SetBins(80,80);//xsec-binning is 25 GeV
365     TH2F * plane = Area<SUSY_XSECS>( xsects, total_xsec );
366     SetBins(400,400);//mass-binning is 5 GeV
367     plane->SetTitle("The mSUGRA parameter space;m_{0} [GeV];m_{1/2} [GeV];#sigma_{SUSY incl.} [pb]");
368     plane->GetYaxis()->SetTitleOffset( 1.8 );
369     plane->GetZaxis()->SetTitleOffset( 1.1 );
370     plane->GetZaxis()->SetRangeUser(0.00001, 1000.);
371     gPad->SetLogz(1);
372     plane->Draw("COLZ");
373    
374     ///Now, add what ever should be displayed else
375     //---------------------------------------------------------------------------//
376     TLatex t;
377     t.SetTextColor(14);
378     t.SetTextAlign(12);
379     t.SetTextSize(0.02);
380    
381     ///gluino
382     gStyle->SetLineColor(15);//gray
383     MyTGraph * m_gluino5 = IsoMassLine( points, gluino, 500. ); m_gluino5->Draw();
384     MyTGraph * m_gluino10 = IsoMassLine( points, gluino, 1000. ); m_gluino10->Draw();
385     MyTGraph * m_gluino20 = IsoMassLine( points, gluino, 2000. ); m_gluino20->Draw();
386     MyTGraph * m_gluino30 = IsoMassLine( points, gluino, 3000. ); m_gluino30->Draw();
387     MyTGraph * m_gluino40 = IsoMassLine( points, gluino, 4000. ); m_gluino40->Draw();
388     t.DrawLatex(360.,450.,"#tilde{g}=1000 GeV");
389     t.DrawLatex(600.,930.,"#tilde{g}=2000 GeV");
390     t.DrawLatex(870.,1420.,"#tilde{g}=3000 GeV");
391     t.DrawLatex(1150.,1930.,"#tilde{g}=4000 GeV");
392    
393     ///neutralino
394     //MyTGraph * m_neutralino1 = IsoMassLine( neutralino_1, 500. ); m_neutralino1->Draw();
395    
396     ///left-handed up-like squarks
397     MyTGraph * m_ul5 = IsoMassLine( points, ul, 500. ); m_ul5->Draw();
398     MyTGraph * m_ul10 = IsoMassLine( points, ul, 1000. ); m_ul10->Draw();
399     MyTGraph * m_ul15 = IsoMassLine( points, ul, 1500. ); m_ul15->Draw();
400     MyTGraph * m_ul20 = IsoMassLine( points, ul, 2000. ); m_ul20->Draw();
401     MyTGraph * m_ul25 = IsoMassLine( points, ul, 2500. ); m_ul25->Draw();
402     MyTGraph * m_ul30 = IsoMassLine( points, ul, 3000. ); m_ul30->Draw();
403     MyTGraph * m_ul35 = IsoMassLine( points, ul, 3500. ); m_ul35->Draw();
404     MyTGraph * m_ul40 = IsoMassLine( points, ul, 4000. ); m_ul40->Draw();
405     t.DrawLatex(800.,250.,"#tilde{u}_{L}=1000 GeV");
406     t.DrawLatex(1670.,470.,"#tilde{u}_{L}=2000 GeV");
407     t.DrawLatex(1670.,1250.,"#tilde{u}_{L}=3000 GeV");
408     t.DrawLatex(1670.,1930.,"#tilde{u}_{L}=4000 GeV");
409    
410     ///Contours
411     /*
412     if (contour_no_solution) {
413     //contour_no_solution->Add(max_x+5,430);
414     //contour_no_solution->Add(max_x+5,min_y-5);
415     //contour_no_solution->Add(min_x,min_y-5);
416     contour_no_solution->SetFillColor( 1 );
417     contour_no_solution->SetFillStyle( 3454 );
418     contour_no_solution->Draw("lf");
419     contour_no_solution->Draw();
420     }
421     */
422     if (contour_charged_LSP){
423     //contour_charged_LSP->Add(min_x-5, 726.25);
424     //contour_charged_LSP->Add(min_x-5, max_y+5);
425     contour_charged_LSP->SetFillColor( 6 );
426     contour_charged_LSP->SetFillStyle( 3454 );
427     contour_charged_LSP->SetLineColor( 6 );
428     contour_charged_LSP->SetLineWidth( 4 );
429     contour_charged_LSP->Draw("lf");
430     contour_charged_LSP->Draw();
431     }
432     t.SetTextSize(0.03);t.SetTextColor( 6 );
433     t.DrawLatex(200.,1800.,"charged LSP (#tilde{#tau})");
434    
435     // m(gluino)==m(squark)
436     gStyle->SetLineColor(1);//black
437     MyTGraph * m_ulgl = IsoMassLine( points, sq_equal_gluino, 0.2 );
438     m_ulgl->SetLineWidth(3); m_ulgl->Draw();
439     t.SetTextSize(0.03);t.SetTextColor( 1 );
440     t.DrawLatex( 750.,1200.,"m(#tilde{u}_{L}) < m(#tilde{g})");
441     t.DrawLatex(1350.,1150.,"m(#tilde{u}_{L}) > m(#tilde{g})");
442    
443    
444     // xs(gluino)==xs(squark)
445     MyTGraph * xs_ulgl = IsoMassLine( xsects, xs_sq_equal_gluino, 5.0 );
446     xs_ulgl->SetLineWidth(3); xs_ulgl->Draw();
447     t.SetTextSize(0.02);t.SetTextColor( 1 );
448     t.DrawLatex( 890.,510.,"#sigma(#tilde{q}#tilde{q}) > #sigmam(#tilde{g}#tilde{g})");
449     t.DrawLatex(1280.,450.,"#sigma(#tilde{q}#tilde{q}) < #sigmam(#tilde{g}#tilde{g})");
450    
451     TMarker m;
452     m.SetMarkerStyle(29);
453     m.SetMarkerColor(2);
454     m.SetMarkerSize(1.5);
455     //m.DrawMarker(1000,1000);
456     //DONE
457     c1->SaveAs("dirt.eps");
458     gApplication->Run(); //interactive
459    
460     //ps.NewPage();
461    
462     }
463    
464     int scan(void)
465     {
466     gStyle->SetHistFillColor(0);
467     gStyle->SetPalette(1);
468     //gStyle->SetFillColor(0);
469     //gStyle->SetOptStat(kFALSE);
470     gStyle->SetCanvasColor(0);
471     gStyle->SetCanvasBorderMode(0);
472     gStyle->SetPadColor(0);
473     gStyle->SetPadBorderMode(0);
474     gStyle->SetFrameBorderMode(0);
475    
476     gStyle->SetTitleFillColor(0);
477     gStyle->SetTitleBorderSize(0);
478     gStyle->SetTitleX(0.10);
479     gStyle->SetTitleY(0.98);
480     gStyle->SetTitleW(0.8);
481     gStyle->SetTitleH(0.06);
482    
483     gStyle->SetErrorX(0);
484     gStyle->SetStatColor(0);
485     gStyle->SetStatBorderSize(0);
486     gStyle->SetStatX(0);
487     gStyle->SetStatY(0);
488     gStyle->SetStatW(0);
489     gStyle->SetStatH(0);
490    
491     gStyle->SetTitleFont(22);
492     gStyle->SetLabelFont(22,"X");
493     gStyle->SetLabelFont(22,"Y");
494     gStyle->SetLabelFont(22,"Z");
495     gStyle->SetLabelSize(0.03,"X");
496     gStyle->SetLabelSize(0.03,"Y");
497     gStyle->SetLabelSize(0.03,"Z");
498    
499     //gROOT->SetStyle("MyStyle");
500    
501    
502     TScan * scan = new TScan();
503     scan->DoStuff();
504     delete scan;
505     }
506    
507     int main(int argc, char** argv)
508     {
509     TRint *theApp = new TRint("ROOT example", &argc, argv);
510     //theApp->Run();
511     //TApplication theApp("App", &argc, argv);
512    
513     return scan();
514     }