CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
bin_peak2D.cxx
Go to the documentation of this file.
1#include <TChain.h>
2#include <TFile.h>
3#include <vector>
4#include <math.h>
5#include "TH1.h"
6#include "TH2.h"
7#include "TGraph.h"
8#include "TCanvas.h"
9#include <iostream>
10#include <fstream>
11#include <string>
12#include <sstream>
13#include <typeinfo.h>
14using namespace std;
15
16
17
18int bin_peak2D()
19{
20 double M_PI = 3.1415926;
21//////////////////////////////////////////////////////////
22// This file has been automatically generated
23// (Wed Jan 10 14:20:10 2018 by ROOT version5.34/09)
24// from TTree hit/hit
25// found on file: hough_hist_test.root
26//////////////////////////////////////////////////////////
27
28
29//Reset ROOT and connect tree file
30 //gROOT->Reset();
31 //TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("hough_hist_test.root");
32 //if (!f) {
33 //f = new TFile("hough_hist_test.root");
34 //}
35 //TTree* tree;
36 //f->GetObject("hit",tree);
37
38 TChain * hit = new TChain("hit");
39 hit->Add("hough_hist_test.root");
40//Declaration of leaves types
41 Int_t hit_run;
42 Int_t hit_evt;
43 Int_t hit_nhit;
44 Int_t hit_hitid[1000];
45 Int_t hit_layer[1000];
46 Int_t hit_wire[1000];
47 Double_t hit_x[1000];
48 Double_t hit_y[1000];
49 Double_t hit_z[1000];
50 Double_t hit_driftdist[1000];
51 Double_t hit_drifttime[10000];
52 Int_t hit_flag[1000];
53 Double_t hit_truth_x[1000];
54 Double_t hit_truth_y[1000];
55 Double_t hit_truth_z[1000];
56 Double_t hit_truth_drift[1000];
57 Int_t hit_truth_ambig[1000];
58
59 // Set branch addresses.
60 hit->SetBranchAddress("hit_run",&hit_run);
61 hit->SetBranchAddress("hit_evt",&hit_evt);
62 hit->SetBranchAddress("hit_nhit",&hit_nhit);
63 hit->SetBranchAddress("hit_hitid",hit_hitid);
64 hit->SetBranchAddress("hit_layer",hit_layer);
65 hit->SetBranchAddress("hit_wire",hit_wire);
66 hit->SetBranchAddress("hit_x",hit_x);
67 hit->SetBranchAddress("hit_y",hit_y);
68 hit->SetBranchAddress("hit_z",hit_z);
69 hit->SetBranchAddress("hit_driftdist",hit_driftdist);
70 hit->SetBranchAddress("hit_drifttime",hit_drifttime);
71 hit->SetBranchAddress("hit_flag",hit_flag);
72 hit->SetBranchAddress("hit_truth_x",hit_truth_x);
73 hit->SetBranchAddress("hit_truth_y",hit_truth_y);
74 hit->SetBranchAddress("hit_truth_z",hit_truth_z);
75 hit->SetBranchAddress("hit_truth_drift",hit_truth_drift);
76 hit->SetBranchAddress("hit_truth_ambig",hit_truth_ambig);
77
78// This is the loop skeleton
79// To read only selected branches, Insert statements like:
80// hit->SetBranchStatus("*",0); // disable all branches
81// TTreePlayer->SetBranchStatus("branchname",1); // activate branchname
82
83 TAxis *x_axis,*y_axis;
84 stringstream ssname;
85 string sname;
86 const char * name;
87 double peak_rate_mean = 0;
88 double peak_height_rate_mean = 0;
89 double peak_rate_RMS = 0;
90 double peak_height_rate_RMS = 0;
91 TGraphErrors *gr_peak_rate_mean = new TGraphErrors();
92 TGraphErrors *gr_peak_height_rate_mean = new TGraphErrors();
93 int ibin = 0;
94 int nbin = 100 ;
95 while(nbin <= 3000){
96 int x_bin = nbin;
97 int y_bin = nbin;
98 double x_max = M_PI;
99 double y_max = 0.1;
100 TH1D *h_peak_rate = new TH1D("peak_rate","peak_rate", 40,0,2.0);
101 TH1D *h_peak_height_rate = new TH1D("peak_height_rate","peak_height_rate", 40,0,2.0);
102 TH1D *h_layersOnPeak = new TH1D("layersOnPeak","layersOnPeak",50,0,50);
103 TH1D *h_layersOfEvent = new TH1D("layersOfEvent","layersOfEvent",50,0,50);
104
105 Long64_t nentries = hit->GetEntries();
106 int event =11;
107 Long64_t nbytes = 0;
108 for (Long64_t i=0; i<nentries;i++)
109 //for (Long64_t i=0; i<1000;i++)
110 //for (Long64_t i=event; i<event+1;i++)
111 {
112 nbytes += hit->GetEntry(i);
113 TH2D* hough = new TH2D("hough","hough",x_bin,0,x_max,y_bin,-y_max,y_max);
114 TH2D *houghhit[100];
115 int layersOnPeak[100];
116 int nhit(0),ihit(0);
117 for(int j=0;j<hit_nhit;j++){
118 if(hit_flag[j] !=0)continue;
119 nhit++;
120 //if(hit_layer[j] >35)continue;
121 //cout<<hit_layer[j]<<endl;
122 //fill_histogram(hit_x[j],hit_y[j],hit_driftdist[j],x_bin*2,hough);
123 ssname.str("");
124 ssname <<"layer:"<<hit_layer[j]<<" ,wire"<<hit_wire[j];
125 sname = ssname.str();
126 name = sname.c_str();
127 layersOnPeak[ihit] = hit_layer[j];
128 houghhit[ihit] = new TH2D(name,name,x_bin,0,x_max,y_bin,-y_max,y_max);
129 fill_histogram(hit_x[j],hit_y[j],hit_driftdist[j],x_bin*2,houghhit[ihit]);
130 hough->Add(houghhit[ihit]);
131 ihit++;
132 h_layersOfEvent->Fill(hit_layer[j]);
133 }
134 int binx,biny,binz;
135 hough->GetMaximumBin(binx,biny,binz);
136 int height = hough->GetBinContent(binx,biny);
137 //double theta = hough->GetXaxis()->GetBinCenter(binx);
138 //double rho = hough->GetYaxis()->GetBinCenter(biny);
139
140 ///*
141 int nhot(0);
142 //TCanvas *c = new TCanvas("houghspace","houghspace",1000,1000 );
143 for(int ihit=0;ihit<nhit;ihit++){
144 //ssname.str("");
145 //ssname <<ihit<<".pdf";
146 //sname = ssname.str();
147 //name = sname.c_str();
148 //c = new TCanvas(name,name,700,500);
149 //houghhit[ihit]->Draw();
150 //cout<<houghhit[ihit]->GetBinContent(binx,biny)<<endl;
151 if(houghhit[ihit]->GetBinContent(binx,biny)>0){
152 nhot++;
153 h_layersOnPeak->Fill(layersOnPeak[ihit]);
154 }
155 //cout<<layersOnPeak[ihit]<<" "<<houghhit[ihit]->GetBinContent(binx,biny)<<endl;
156 //c->SaveAs(name);
157 delete houghhit[ihit];
158 }
159 //cout<<nhit<<" "<<height<<" "<<nhot<<" "<<(double)height/(double)nhit<<" "<<(double)nhot/(double)nhit<<endl;
160 double peak_rate = (double)nhot/(double)nhit;
161 double peak_height_rate = (double)height/(double)nhit;
162 h_peak_rate->Fill(peak_rate);
163 h_peak_height_rate->Fill(peak_height_rate);
164
165 //TCanvas *c1 = new TCanvas("houghspace","houghspace",1000,500 );
166 //c1->Divide(2,1);
167 //c1->cd(1);
168 //hough->Draw();
169 //c1->cd(2);
170 //hough->Draw("LEGO2Z");
171 //c1->SaveAs("houghspace.pdf");
172
173 delete hough;
174 //c1->SaveAs("houghspace.pdf");
175 //*/
176 }
177 ///*
178 ssname.str("");
179 ssname <<x_bin<<"bins_layersOnPeak";
180 sname = ssname.str();
181 name = sname.c_str();
182 TCanvas *c2 = new TCanvas(name,name,700,500 );
183 h_layersOnPeak->Draw("same");
184 h_layersOnPeak->SetLineColor(4);
185 x_axis = h_layersOnPeak->GetXaxis();
186 y_axis = h_layersOnPeak->GetYaxis();
187 x_axis->SetTitle("Layer");
188 x_axis->CenterTitle();
189 y_axis->SetTitle("Entries");
190 y_axis->CenterTitle();
191 yMax = h_layersOnPeak->GetMaximum() > h_layersOfEvent->GetMaximum()? h_layersOnPeak->GetMaximum():h_layersOfEvent->GetMaximum();
192 h_layersOnPeak->GetYaxis()->SetRangeUser(0,yMax*1.2);
193 h_layersOfEvent->Draw("same");
194 h_layersOfEvent->SetLineColor(2);
195 TLegend* leg2= new TLegend(0.12 ,0.80,0.33,0.89);
196 leg2->AddEntry(h_layersOnPeak,"hits on peak","l");
197 leg2->AddEntry(h_layersOfEvent,"hits of event","l");
198 leg2->SetFillColor(kWhite);
199 leg2->SetLineColor(kWhite);
200 leg2->Draw();
201 ssname.str("");
202 ssname <<x_bin<<"bins_layersOnPeak"<<".pdf";
203 sname = ssname.str();
204 name = sname.c_str();
205 c2->SaveAs(name);
206 ssname.str("");
207 ssname <<x_bin<<"bins_layersOnPeak"<<".eps";
208 sname = ssname.str();
209 name = sname.c_str();
210 c2->SaveAs(name);
211 ssname.str("");
212 ssname <<x_bin<<"bins_layersOnPeak"<<".png";
213 sname = ssname.str();
214 name = sname.c_str();
215 c2->SaveAs(name);
216
217 peak_rate_mean = h_peak_rate->GetMean();
218 peak_height_rate_mean = h_peak_height_rate->GetMean();
219 peak_rate_RMS = h_peak_rate->GetRMS();
220 peak_height_rate_RMS = h_peak_height_rate->GetRMS();
221 cout<<x_bin<<" "<<peak_rate_mean<<" "<<peak_rate_RMS<<" "<<peak_height_rate_mean<<" "<<peak_height_rate_RMS<<endl;
222 gr_peak_rate_mean->SetPoint(ibin,x_bin,peak_rate_mean);
223 gr_peak_rate_mean->SetPointError(ibin,0,peak_rate_RMS);
224 gr_peak_height_rate_mean->SetPoint(ibin,x_bin,peak_height_rate_mean);
225 gr_peak_height_rate_mean->SetPointError(ibin,0,peak_height_rate_RMS);
226 ssname.str("");
227 ssname <<x_bin<<"bins_rate";
228 sname = ssname.str();
229 name = sname.c_str();
230 TCanvas *c3 = new TCanvas(name,name,700,500 );
231 h_peak_rate->Draw();
232 h_peak_height_rate->Draw("same");
233 h_peak_rate->SetLineColor(4);
234 h_peak_height_rate->SetLineColor(2);
235 x_axis = h_peak_rate->GetXaxis();
236 y_axis = h_peak_rate->GetYaxis();
237 //x_axis->SetLimits(0,80 );
238 x_axis->SetTitle("Rate");
239 x_axis->CenterTitle();
240 //x_axis->SetTitleSize(0.05);
241 //x_axis->SetLabelSize(0.04);
242 //x_axis->SetTitleOffset(1.05);
243 y_axis->SetTitle("Entries");
244 y_axis->CenterTitle();
245 //y_axis->SetTitleSize(0.05);
246 //y_axis->SetLabelSize(0.04);
247 //y_axis->SetTitleOffset(1.25);
248 yMax = h_peak_rate->GetMaximum() > h_peak_height_rate->GetMaximum()? h_peak_rate->GetMaximum():h_peak_height_rate->GetMaximum();
249 h_peak_rate->GetYaxis()->SetRangeUser(0,yMax*1.2);
250 TLegend* leg3= new TLegend(0.12 ,0.80,0.33,0.89);
251 leg3->AddEntry(h_peak_rate,"hit on peak","l");
252 leg3->AddEntry(h_peak_height_rate,"peak height","l");
253 leg3->SetFillColor(kWhite);
254 leg3->SetLineColor(kWhite);
255 leg3->Draw();
256 ssname.str("");
257 ssname <<x_bin<<"bins_rate"<<".pdf";
258 sname = ssname.str();
259 name = sname.c_str();
260 c3->SaveAs(name);
261 ssname.str("");
262 ssname <<x_bin<<"bins_rate"<<".eps";
263 sname = ssname.str();
264 name = sname.c_str();
265 c3->SaveAs(name);
266 ssname.str("");
267 ssname <<x_bin<<"bins_rate"<<".png";
268 sname = ssname.str();
269 name = sname.c_str();
270 c3->SaveAs(name);
271 delete h_layersOnPeak;
272 delete h_peak_rate;
273 delete h_peak_height_rate;
274
275 int bin;
276 if(nbin<500) bin = 50;
277 else if(nbin<1000) bin =100;
278 else if(nbin<2000)bin =200;
279 else if(nbin<5000)bin =500;
280 else bin = 1000;
281 nbin += bin;
282 ibin++;
283 //*/
284 }
285 ///*
286 TCanvas *c4 = new TCanvas("peak_rate_mean","peak_rate_mean",700,500 );
287 gr_peak_rate_mean->Draw("AP");
288 x_axis = gr_peak_rate_mean->GetXaxis();
289 y_axis = gr_peak_rate_mean->GetYaxis();
290 //x_axis->SetLimits(0,80 );
291 x_axis->SetTitle("nbin");
292 x_axis->CenterTitle();
293 //x_axis->SetTitleSize(0.05);
294 //x_axis->SetLabelSize(0.04);
295 //x_axis->SetTitleOffset(1.05);
296 y_axis->SetTitle("mean rate");
297 y_axis->CenterTitle();
298 //y_axis->SetTitleSize(0.05);
299 //y_axis->SetLabelSize(0.04);
300 //y_axis->SetTitleOffset(1.25);
301 gr_peak_rate_mean->SetMinimum(0.0);
302 gr_peak_rate_mean->SetMaximum(2.0);
303 gr_peak_rate_mean->SetMarkerSize(1.0);
304 gr_peak_rate_mean->SetMarkerStyle(24);
305 gr_peak_rate_mean->SetMarkerColor(4);
306 gr_peak_height_rate_mean->Draw("Psame");
307 gr_peak_height_rate_mean->SetMarkerSize(1.0);
308 gr_peak_height_rate_mean->SetMarkerStyle(24);
309 gr_peak_height_rate_mean->SetMarkerColor(2);
310 TLegend* leg4= new TLegend(0.60 ,0.75,0.85,0.85);
311 leg4->SetFillColor(kWhite);
312 leg4->SetLineColor(kWhite);
313 leg4->AddEntry(gr_peak_rate_mean,"hit on peak","p");
314 leg4->AddEntry(gr_peak_height_rate_mean,"peak height","p");
315 leg4->Draw();
316
317 c4->SaveAs("peak_rate_mean.eps");
318 c4->SaveAs("peak_rate_mean.png");
319 c4->SaveAs("peak_rate_mean.pdf");
320 //*/
321}
322
323
324void fill_histogram(double x, double y, double r, int nbin, TH2D* hough)
325{
326 double M_PI = 3.1415926;
327 double U = 2*x/(x*x+y*y-r*r);
328 double V = 2*y/(x*x+y*y-r*r);
329 double R = 2*r/(x*x+y*y-r*r);
330 double delta_phi = 2.*M_PI/nbin;
331 double phi = -delta_phi/2;
332 for(int i =0; i<nbin; i++){
333 phi += delta_phi;
334 double u = U + R*cos(phi);
335 double v = V + R*sin(phi);
336 double normal = (v-V)/(u-U);
337 double k = -1./normal;
338 double b = v - k*u;
339 double x_cross = -b/(k+1/k);
340 double y_cross = b/(1+k*k);
341 // std::cout<<"x y centerX centerY "<<x<<" "<<y<<" "<<centerX<<" "<<centerY<<std::endl;
342 // std::cout<<"k b "<<k<<" "<<b<<std::endl;
343 // std::cout<<"xcross ycross "<<x_cross<<" "<<y_cross<<std::endl;
344
345 double rho=sqrt(x_cross*x_cross+y_cross*y_cross);
346 double theta=atan2(y_cross,x_cross);
347 if(theta<0) {
348 theta=theta+M_PI;
349 rho=-rho;
350 }
351 if( normal ==0 && x>0) {
352 rho = fabs(x);
353 theta = 0;
354 }
355 if( normal ==0 && x<0) {
356 rho = -fabs(x);
357 theta = M_PI;
358 }
359 hough->Fill(theta,rho);
360 }
361 //double slant = _y*cos(_theta)-_x*sin(_theta);
362 //_slant = slant;
363// std::cout<<"THETA RHO "<<theta_temp<<" "<<rho_temp<<std::endl;
364// std::cout<<std::endl;
365}
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
nentries
Int_t nentries
Definition: DataBase/tau_mode.c:8
bin
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per bin
Definition: FoamA.h:85
sin
double sin(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:213
v
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
Definition: KarLud.h:35
M_PI
#define M_PI
Definition: TConstant.h:4
bin_peak2D
int bin_peak2D()
Definition: bin_peak2D.cxx:18
fill_histogram
void fill_histogram(double x, double y, double r, int nbin, TH2D *hough)
Definition: bin_peak2D.cxx:324