1 ///
2 /// Use this file from root. First compile it and then run showjets(...)
3 /// on a file produced by ClusterSequence::print_jets_for_root(...)
4 ///
5 /// root> .L jet-plots.C+
6 /// root> showjets(filename [, label])
7 ///
8 #include<vector>
9 #include<fstream>
10 #include<iostream>
11 #include<sstream>
12 #include<cstdlib>
13 #include "TLatex.h"
14
15 #include "TROOT.h"
16 #include "TCanvas.h"
17 #include "TH2.h"
18 #include "THStack.h"
19 #include "TStyle.h"
20 #include "TPaveLabel.h"
21 #include "TColor.h"
22
23 using namespace std;
24
25
26 class JetHist {
27 private:
28 vector<TH2D *> _jets;
29 TH2D * _background;
30 string _comment;
31 public:
32 static double default_etamax;
33 static int default_nbins;
34
35 JetHist(const string & filename, double etamax=default_etamax, int nbins=default_nbins);
36 ~JetHist();
comment()37 string comment() {return _comment;}
38 THStack stack;
jet(int i)39 TH2D * jet(int i) {return i>= 0 ? _jets[i] : _background;}
40
41 };
42
43 double JetHist::default_etamax = 6.0;
44 int JetHist::default_nbins = 40; // y: 2*nbins; phi: nbins
45
46 // get jet "histograms" from filename which is expected to be made of repeated
47 // blocks as follows:
48 // jet# eta phi pt ...
49 // ipart eta phi pt
50 // ipart eta phi pt
51 // ...
52 // #END
JetHist(const string & filename,double etamax,int nbins)53 JetHist::JetHist (const string & filename, double etamax, int nbins) {
54 ifstream file(filename.c_str());
55 string line;
56 //double etamax=6;
57 //double etamax=5;
58 double phimax = 2*3.14159265;
59 //int nbins=40;
60
61 // construct a histogram for the background to the jets
62 ostringstream bname;
63 bname << filename <<"-background";
64 _background = new TH2D(bname.str().c_str(),bname.str().c_str(),
65 2*nbins,-etamax,etamax,nbins,0.0,phimax);
66 //_background = new TH2D(bname.str().c_str(),bname.str().c_str(),
67 // 2*nbins,-etamax,etamax,2,0.0,phimax);
68 _background->SetFillColor(kWhite);
69 // these were supposed to have labelled the axes, but it doesn't work.
70 _background->GetXaxis()->SetTitle("#eta");
71 _background->GetYaxis()->SetTitle("#phi");
72 _background->GetZaxis()->SetTitle("p_{#perp}");
73 stack.Add(_background);
74
75 while (getline(file,line)) {
76 if (line.substr(0,1) != " ") { // all interesting lines start with space?
77 // extract a comment if there is one
78 if (line.substr(0,2) == "# ") {
79 _comment = line.substr(2,line.length()-2);
80 }
81 continue;
82 }
83 ostringstream name;
84 name << filename<<"-jet-"<< _jets.size();
85 TH2D * hist = new TH2D(name.str().c_str(),name.str().c_str(),
86 2*nbins,-etamax,etamax,nbins,0.0,phimax);
87 int i;
88 double eta, phi, pt;
89 //cout << filename <<": jet "<<_jets.size()<<endl;
90 bool have_line = true;
91 while (have_line || getline(file,line)) {
92 have_line = false;
93 if (line.substr(0,4) == "#END") {break;}
94 istringstream sline(line);
95 sline >> i >> eta >> phi >> pt;
96 //cout << i << " "<<eta<<" "<<phi<<" "<<pt<<endl;
97 hist->Fill(eta,phi,pt); // fill at phi,eta with weight pt
98
99 // workaround for bug in stacks: fill all lower elements of the stack
100 // with a fake amount -- this, miraculously will lead to correct coloring
101 // of the top of the stack!
102 //for (unsigned int j = 0; j < _jets.size(); j++) {
103 // _jets[j]->Fill(phi,eta,1e-7);
104 //}
105 }
106 // give it a colour (whatever that means...)
107 //hist->SetFillColor(_jets.size());
108 int njet = _jets.size();
109 //hist->SetFillColor(njet+2);
110 hist->SetFillColor(njet%50+2); // %50 seems tomake to diff to many-jet case
111 //if (njet == 0) hist->SetFillColor(kRed);
112 //else if (njet == 1) hist->SetFillColor(kBlue);
113 //else hist->SetFillColor(kGreen);
114
115 // add it to the list of jets (so we can delete it later...)
116 _jets.push_back(hist);
117
118 // put it onto the stack
119 stack.Add(hist);
120 //if (njet == 2) break;
121 //break;
122 }
123 }
124
125 // clean up --------------
~JetHist()126 JetHist::~JetHist () {
127 for (unsigned int i = 0; i < _jets.size(); i++) {
128 delete _jets[i];
129 }
130 delete _background;
131 }
132
133 //----------------------------------------------------------------------
134 /// set up a reasonable bunch of colours
set_default_colours(TCanvas * lego)135 void set_default_colours(TCanvas * lego) {
136 gStyle->SetOptStat(0);
137 gStyle->SetPalette(1);
138 gStyle->SetCanvasColor(1);
139 gStyle->SetFrameFillColor(0);
140
141 Int_t cancolor = 0;
142 lego->SetFillColor(cancolor);
143
144 int ngrey = 3;
145 for (int ir = 0; ir < ngrey; ir++) {
146 for (int ig = 0; ig < ngrey; ig++) {
147 for (int ib = 0; ib < ngrey; ib++) {
148 int icol = 7+ir + ngrey *ig + ngrey*ngrey * ib;
149 TColor * color=(TColor*)(gROOT->GetListOfColors()->At(icol));
150 if (icol == 7) {
151 // avoid white -- put grey instead
152 color->SetRGB(0.5,0.5,0.5);
153 } else {
154 color->SetRGB(1-ir*1.0/ngrey,1-ig*1.0/ngrey,1-ib*1.0/ngrey);
155 }
156 }
157 }
158 }
159 }
160
161 //----------------------------------------------------------------------
162 /// show the jets contained in filename (as produced by
163 /// ClusterSequence::print_jets_for_root()), with an optional label
164 TCanvas * showjets (const char * filename, const char * label = 0) {
165
166 // display the various 2-d drawing options
167 gROOT->Reset();
168
169 // set up canvas
170 TCanvas * lego = new TCanvas("lego","lego options",400,50,800,600);
171 lego->SetTheta(30.0);
172 lego->SetPhi(20.0);
173
174 // orientation used for plots in subtraction paper
175 //lego->SetTheta(62.15);
176 //lego->SetPhi(9.15);
177
178 ////vector<double> col
179
180 set_default_colours(lego);
181
182 TPaveLabel pl;
183
184 JetHist * jets = new JetHist(filename);
185 jets->stack.Draw("lego1"); // cyl does not work with 5.16
186 if (label != 0) {
187 Float_t x1=0.63, y1=0.875, x2=0.95, y2=0.925;
188 pl.DrawPaveLabel(x1,y1,x2,y2,label,"brNDC");
189 } else if (jets->comment() != "") {
190 Float_t x1=0.15, y1=0.875, x2=0.95, y2=0.925;
191 pl.DrawPaveLabel(x1,y1,x2,y2,jets->comment().c_str(),"brNDC");
192 }
193
194 // normal histogram labels not working, so draw them by hand
195 TLatex l;
196 l.SetTextAlign(22);
197 l.SetTextSize(0.05);
198 //l.DrawLatex(0.0,0.85,"anti-k_{t}, R=1");
199
200 l.SetTextSize(0.04);
201 l.DrawLatex(0.20,-0.98,"y");
202 l.SetTextAlign(32);
203 l.DrawLatex(-0.7,0.8,"p_{t} [GeV]");
204 l.DrawLatex(-0.6,-0.78,"#phi");
205
206 // do not delete jets -- otherwise you lose everything!;
207
208 return lego;
209 ///
210 lego->Update();
211 }
212