CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
compare_trk_par.cxx
Go to the documentation of this file.
1 /*
2 this script is used to compare the reconstrcuted track parameters between MdcHoughFinder and ROOT
3 */
4int compare_trk_par()
5{
6 //////////////////////////////////////////////////////////
7 // This file has been automatically generated
8 // (Wed Jan 10 14:20:51 2018 by ROOT version5.34/09)
9 // from TTree trk/trk
10 // found on file: hough_hist_test.root
11 //////////////////////////////////////////////////////////
12
13
14 //Reset ROOT and connect tree file
15 //gROOT->Reset();
16 //TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("hough_hist_test.root");
17 //if (!f) {
18 //f = new TFile("hough_hist_test.root");
19 //}
20 //TTree* tree;
21 //f->GetObject("trk",tree);
22 const double M_PI = 3.1415926;
23
24 TChain * trk = new TChain("trk");
25 trk->Add("hough_hist_test.root");
26 //Declaration of leaves types
27 Int_t trk_run;
28 Int_t trk_evt;
29 Int_t trk_ntrk;
30 Int_t trk_trackId[100];
31 Int_t trk_nhop[100];
32 Int_t trk_nhot[100];
33 Double_t trk_dr[100];
34 Double_t trk_phi0[100];
35 Double_t trk_kappa[100];
36 Double_t trk_dz[100];
37 Double_t trk_tanl[100];
38 Double_t trk_Xc[100];
39 Double_t trk_Yc[100];
40 Double_t trk_R[100];
41 Double_t trk_p[100];
42 Double_t trk_pt[100];
43 Double_t trk_px[100];
44 Double_t trk_py[100];
45 Double_t trk_pz[100];
46 Int_t trk_q[100];
47 Double_t trk_truth_dr[100];
48 Double_t trk_truth_phi0[100];
49 Double_t trk_truth_kappa[100];
50 Double_t trk_truth_dz[100];
51 Double_t trk_truth_tanl[100];
52 Double_t trk_truth_Xc[100];
53 Double_t trk_truth_Yc[100];
54 Double_t trk_truth_R[100];
55 Double_t trk_truth_cosTheta[100];
56 Double_t trk_truth_p[100];
57 Double_t trk_truth_pt[100];
58 Double_t trk_truth_px[100];
59 Double_t trk_truth_py[100];
60 Double_t trk_truth_pz[100];
61 Int_t trk_truth_q[100];
62
63 // Set branch addresses.
64 trk->SetBranchAddress("trk_run",&trk_run);
65 trk->SetBranchAddress("trk_evt",&trk_evt);
66 trk->SetBranchAddress("trk_ntrk",&trk_ntrk);
67 trk->SetBranchAddress("trk_trackId",trk_trackId);
68 trk->SetBranchAddress("trk_nhop",trk_nhop);
69 trk->SetBranchAddress("trk_nhot",trk_nhot);
70 trk->SetBranchAddress("trk_dr",trk_dr);
71 trk->SetBranchAddress("trk_phi0",trk_phi0);
72 trk->SetBranchAddress("trk_kappa",trk_kappa);
73 trk->SetBranchAddress("trk_dz",trk_dz);
74 trk->SetBranchAddress("trk_tanl",trk_tanl);
75 trk->SetBranchAddress("trk_Xc",trk_Xc);
76 trk->SetBranchAddress("trk_Yc",trk_Yc);
77 trk->SetBranchAddress("trk_R",trk_R);
78 trk->SetBranchAddress("trk_p",trk_p);
79 trk->SetBranchAddress("trk_pt",trk_pt);
80 trk->SetBranchAddress("trk_px",trk_px);
81 trk->SetBranchAddress("trk_py",trk_py);
82 trk->SetBranchAddress("trk_pz",trk_pz);
83 trk->SetBranchAddress("trk_q",trk_q);
84 trk->SetBranchAddress("trk_truth_dr",trk_truth_dr);
85 trk->SetBranchAddress("trk_truth_phi0",trk_truth_phi0);
86 trk->SetBranchAddress("trk_truth_kappa",trk_truth_kappa);
87 trk->SetBranchAddress("trk_truth_dz",trk_truth_dz);
88 trk->SetBranchAddress("trk_truth_tanl",trk_truth_tanl);
89 trk->SetBranchAddress("trk_truth_Xc",trk_truth_Xc);
90 trk->SetBranchAddress("trk_truth_Yc",trk_truth_Yc);
91 trk->SetBranchAddress("trk_truth_R",trk_truth_R);
92 trk->SetBranchAddress("trk_truth_cosTheta",trk_truth_cosTheta);
93 trk->SetBranchAddress("trk_truth_p",trk_truth_p);
94 trk->SetBranchAddress("trk_truth_pt",trk_truth_pt);
95 trk->SetBranchAddress("trk_truth_px",trk_truth_px);
96 trk->SetBranchAddress("trk_truth_py",trk_truth_py);
97 trk->SetBranchAddress("trk_truth_pz",trk_truth_pz);
98 trk->SetBranchAddress("trk_truth_q",trk_truth_q);
99
100 // This is the loop skeleton
101 // To read only selected branches, Insert statements like:
102 // trk->SetBranchStatus("*",0); // disable all branches
103 // TTreePlayer->SetBranchStatus("branchname",1); // activate branchname
104
105 //////////////////////////////////////////////////////////
106 // This file has been automatically generated
107 // (Wed Jan 10 14:20:10 2018 by ROOT version5.34/09)
108 // from TTree hit/hit
109 // found on file: hough_hist_test.root
110 //////////////////////////////////////////////////////////
111
112
113 //Reset ROOT and connect tree file
114 //gROOT->Reset();
115 //TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("hough_hist_test.root");
116 //if (!f) {
117 //f = new TFile("hough_hist_test.root");
118 //}
119 //TTree* tree;
120 //f->GetObject("hit",tree);
121
122 gStyle->SetOptFit(0111);
123 TChain * hit = new TChain("hit");
124 hit->Add("hough_hist_test.root");
125 //Declaration of leaves types
126 Int_t hit_run;
127 Int_t hit_evt;
128 Int_t hit_nhit;
129 Int_t hit_hitid[1000];
130 Int_t hit_layer[1000];
131 Int_t hit_wire[1000];
132 Double_t hit_x[1000];
133 Double_t hit_y[1000];
134 Double_t hit_z[1000];
135 Double_t hit_driftdist[1000];
136 Double_t hit_drifttime[10000];
137 Int_t hit_flag[1000];
138 Double_t hit_truth_x[1000];
139 Double_t hit_truth_y[1000];
140 Double_t hit_truth_z[1000];
141 Double_t hit_truth_drift[1000];
142 Int_t hit_truth_ambig[1000];
143
144 // Set branch addresses.
145 hit->SetBranchAddress("hit_run",&hit_run);
146 hit->SetBranchAddress("hit_evt",&hit_evt);
147 hit->SetBranchAddress("hit_nhit",&hit_nhit);
148 hit->SetBranchAddress("hit_hitid",hit_hitid);
149 hit->SetBranchAddress("hit_layer",hit_layer);
150 hit->SetBranchAddress("hit_wire",hit_wire);
151 hit->SetBranchAddress("hit_x",hit_x);
152 hit->SetBranchAddress("hit_y",hit_y);
153 hit->SetBranchAddress("hit_z",hit_z);
154 hit->SetBranchAddress("hit_driftdist",hit_driftdist);
155 hit->SetBranchAddress("hit_drifttime",hit_drifttime);
156 hit->SetBranchAddress("hit_flag",hit_flag);
157 hit->SetBranchAddress("hit_truth_x",hit_truth_x);
158 hit->SetBranchAddress("hit_truth_y",hit_truth_y);
159 hit->SetBranchAddress("hit_truth_z",hit_truth_z);
160 hit->SetBranchAddress("hit_truth_drift",hit_truth_drift);
161 hit->SetBranchAddress("hit_truth_ambig",hit_truth_ambig);
162
163 // This is the loop skeleton
164 // To read only selected branches, Insert statements like:
165 // hit->SetBranchStatus("*",0); // disable all branches
166 // TTreePlayer->SetBranchStatus("branchname",1); // activate branchname
167
168 //////////////////////////////////////////////////////////
169 // This file has been automatically generated
170 // (Wed Jan 10 14:20:35 2018 by ROOT version5.34/09)
171 // from TTree hot/hot
172 // found on file: hough_hist_test.root
173 //////////////////////////////////////////////////////////
174
175
176 //Reset ROOT and connect tree file
177 //gROOT->Reset();
178 //TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("hough_hist_test.root");
179 //if (!f) {
180 //f = new TFile("hough_hist_test.root");
181 //}
182 //TTree* tree;
183 //f->GetObject("hot",tree);
184
185 TFile* f = TFile::Open("hough_hist_test.root");
186 TChain * hot = new TChain("hot");
187 hot->Add("hough_hist_test.root");
188 //Declaration of leaves types
189 Int_t hot_run;
190 Int_t hot_evt;
191 Int_t hot_trk;
192 Int_t hot_nhot;
193 Int_t hot_hitid[1000];
194 Int_t hot_layer[1000];
195 Int_t hot_wire[1000];
196 Double_t hot_x[1000];
197 Double_t hot_y[1000];
198 Double_t hot_z[1000];
199 Double_t hot_drift[1000];
200 Int_t hot_flag[1000];
201 Double_t hot_deltaD[1000];
202 Double_t hot_truth_x[1000];
203 Double_t hot_truth_y[1000];
204 Double_t hot_truth_z[1000];
205 Double_t hot_truth_drift[1000];
206 Int_t hot_truth_ambig[1000];
207
208 // Set branch addresses.
209 hot->SetBranchAddress("hot_run",&hot_run);
210 hot->SetBranchAddress("hot_evt",&hot_evt);
211 hot->SetBranchAddress("hot_trk",&hot_trk);
212 hot->SetBranchAddress("hot_nhot",&hot_nhot);
213 hot->SetBranchAddress("hot_hitid",hot_hitid);
214 hot->SetBranchAddress("hot_layer",hot_layer);
215 hot->SetBranchAddress("hot_wire",hot_wire);
216 hot->SetBranchAddress("hot_x",hot_x);
217 hot->SetBranchAddress("hot_y",hot_y);
218 hot->SetBranchAddress("hot_z",hot_z);
219 hot->SetBranchAddress("hot_drift",hot_drift);
220 hot->SetBranchAddress("hot_flag",hot_flag);
221 hot->SetBranchAddress("hot_deltaD",hot_deltaD);
222 hot->SetBranchAddress("hot_truth_x",hot_truth_x);
223 hot->SetBranchAddress("hot_truth_y",hot_truth_y);
224 hot->SetBranchAddress("hot_truth_z",hot_truth_z);
225 hot->SetBranchAddress("hot_truth_drift",hot_truth_drift);
226 hot->SetBranchAddress("hot_truth_ambig",hot_truth_ambig);
227
228 // This is the loop skeleton
229 // To read only selected branches, Insert statements like:
230 // hot->SetBranchStatus("*",0); // disable all branches
231 // TTreePlayer->SetBranchStatus("branchname",1); // activate branchname
232 Long64_t nentries = trk->GetEntries();
233 Long64_t nbytes = 0;
234 int event =0 ;
235 //for (Long64_t i=0; i<nentries;i++)
236 //for (Long64_t i=0; i<100 ;i++)
237 for (Long64_t i=event; i<event+1;i++)
238 {
239 nbytes += trk->GetEntry(i);
240 cout<<trk_evt<<" "<<trk_ntrk<<endl;
241 for(int j=0;j<trk_ntrk;j++){
242 TAxis *x_axis,*y_axis;
243 stringstream ssname;
244 string sname;
245 const char * name;
246 int ibin = 0;
247 int nbin =1000 ;
248 while(nbin <= 1000){
249 int x_bin = nbin;
250 int y_bin = nbin;
251 double x_max = M_PI;
252 double y_max = 0.1;
253
254 //int event =1272;
255 Long64_t nentries = hit->GetEntries();
256 Long64_t nbytes = 0;
257 for (Long64_t ii=0; ii<nentries;ii++)
258 //for (Long64_t i=0; i<1000;i++)
259 //for (Long64_t i=event; i<event+1;i++)
260 {
261 nbytes += hit->GetEntry(ii);
262 if(hit_evt!=trk_evt)continue;
263 TH2D* hough = new TH2D("hough","hough",x_bin,0,x_max,y_bin,-y_max,y_max);
264 TH2D *houghhit[1000];
265 int nhit(0),ihit(0);
266 //cout<<i<<" "<<hit_nhit<<endl;
267 for(int jj=0;jj<hit_nhit;jj++){
268 if(hit_flag[jj] !=0)continue;
269 nhit++;
270 //if(hit_layer[jj] >35)continue;
271 //cout<<hit_layer[jj]<<endl;
272 //fill_histogram(hit_x[jj],hit_y[jj],hit_driftdist[jj],x_bin*2,hough);
273 ssname.str("");
274 ssname <<"layer:"<<hit_layer[jj]<<" ,wire"<<hit_wire[jj];
275 sname = ssname.str();
276 name = sname.c_str();
277 houghhit[ihit] = new TH2D(name,name,x_bin,0,x_max,y_bin,-y_max,y_max);
278 fill_histogram(hit_x[jj],hit_y[jj],hit_driftdist[jj],x_bin*2,houghhit[ihit]);
279 hough->Add(houghhit[ihit]);
280 ihit++;
281 }
282
283 //TCanvas *c1 = new TCanvas("houghspace","houghspace",1000,500 );
284 //c1->Divide(2,1);
285 //c1->cd(1);
286 //hough->Draw();
287 //c1->cd(2);
288 //hough->Draw("LEGO2Z");
289
290 int binx,biny,binz;
291 hough->GetMaximumBin(binx,biny,binz);
292 double theta = hough->GetXaxis()->GetBinCenter(binx);
293 double rho = hough->GetYaxis()->GetBinCenter(biny);
294 double t = trk_phi0[j];
295 double r = trk_kappa[j]/333.567;
296 if(t>M_PI)t=t-M_PI;
297 else r=-r;
298 //cout<<"rec: "<<trk_phi0[j]<<" "<<trk_kappa[j]/333.567<<endl;
299 cout<<" rec: "<<t<<setw(15)<<r<<endl;
300 cout<<" map: "<<theta<<setw(15)<<rho<<endl;
301 double dtheta = fabs(theta-t);
302 double drho = fabs(rho-r);
303 double theta_bin = M_PI/100.;
304 double rho_bin = 0.1/100.;
305 cout<<"diff: "<<dtheta<<setw(15)<<drho<<endl;
306 cout<<"nbin: "<<dtheta/theta_bin<<setw(15)<<drho/rho_bin<<endl;
307 cout<<"rec bin: "<<t/theta_bin<<setw(15)<<r/rho_bin<<endl;
308 cout<<"map bin: "<<theta/theta_bin<<setw(15)<<rho/rho_bin<<endl;
309 cout<<endl;
310
311 //theta = hough->GetXaxis()->GetBinCenter(binx);
312 //rho = hough->GetYaxis()->GetBinCenter(biny);
313 for(int ihit=0;ihit<nhit;ihit++){
314 //cout<<layersOnPeak[ihit]<<" "<<houghhit[ihit]->GetBinContent(binx,biny)<<endl;
315 delete houghhit[ihit];
316 }
317 double Xc = cos(theta)/rho;
318 double Yc = sin(theta)/rho;
319 double Rc = 1./fabs(rho);
320 double d0 = sqrt( Xc*Xc + Yc*Yc )-Rc;
321 double phi0 = atan2(Yc,Xc)+M_PI/2.;
322 double omega = 1/Rc;
323 double dr = -d0;
324 phi0 = phi0-M_PI/2;
325 while(phi0>2*M_PI)phi0-=2*M_PI;
326 while(phi0<0)phi0+=2*M_PI;
327 double kappa = -omega*333.567;
328 //cout<<dr<<" "<<phi0<<" "<<kappa<<endl;
329 cout<<hit_evt<<" "<<Xc<<" "<<Yc<<" "<<Rc<<endl;
330 cout<<trk_evt<<" "<<trk_Xc[j]<<" "<<trk_Yc[j]<<" "<<trk_R[j]<<endl;
331 double deltaD(0);
332 for(int jj=0;jj<hit_nhit;jj++){
333 if(hit_flag[jj] !=0)continue;
334 if(hit_layer[jj] <3){
335 double r_cgem = sqrt(hit_x[jj]*hit_x[jj]+hit_y[jj]*hit_y[jj]);
336 double phi_cross = intersect_cylinder(-1,Xc,Yc,r_cgem);
337 double phih_cluster = atan2(hit_y[jj],hit_x[jj]);
338 double dphi = phih_cluster - phi_cross;
339 if(dphi>M_PI)dphi-=2*M_PI;
340 if(dphi<-1*M_PI)dphi+=2*M_PI;
341 deltaD = r_cgem*dphi;
342 }else{
343 double l1l2 = sqrt((hit_x[jj]-Xc)*(hit_x[jj]-Xc)+(hit_y[jj]-Yc)*(hit_y[jj]-Yc));
344 if(l1l2>Rc)deltaD = l1l2-Rc-hit_driftdist[jj];
345 else deltaD = l1l2-Rc+hit_driftdist[jj];
346 }
347 {
348 Long64_t nentries = hot->GetEntries();
349 Long64_t nbytes = 0;
350 for (Long64_t iii=0; iii<nentries;iii++)
351 //for (Long64_t i=0; i<100 ;i++)
352 //for (Long64_t i=event; i<event+1;i++)
353 {
354 nbytes += hot->GetEntry(iii);
355 if(hot_evt!=hit_evt)continue;
356 for(int jjj=0;jjj<hot_nhot;jjj++){
357 if(hot_flag[jjj] !=0)continue;
358 if(hot_layer[jjj]==hit_layer[jj]&&hot_wire[jjj]==hit_wire[jj]&&hot_hitid[jjj]==hit_hitid[jj]){
359 cout<<setw(5)<<hot_layer[jjj]<<setw(15)<<deltaD<<setw(15)<<hot_deltaD[jjj]<<setw(15)<<deltaD-hot_deltaD[jjj]<<endl;
360 }
361 //else cout<<hit_layer[jj]<<endl;
362 //cout<<hot_layer[jjj]<<" "<<hot_deltaD[jjj]<<endl;
363 }
364 }
365 }
366 }
367
368 delete hough;
369 }
370 int bin(0);
371 if(nbin<500) bin = 50;
372 else if(nbin<1000) bin =100;
373 else if(nbin<2000)bin =100;
374 else if(nbin<5000)bin =500;
375 else bin = 1000;
376 nbin += bin;
377 ibin++;
378 }
379 }
380 cout<<endl;
381 }
382}
383
384
385double intersect_cylinder(int charge, double x_center, double y_center, double r_cylinder)
386{
387 const double M_PI = 3.1415926;
388 double phi;
389 double phi_center = atan2(y_center,x_center);
390 double r_center = sqrt(x_center*x_center+y_center*y_center);
391 if(charge==0)return 9999;
392 while(phi_center<0) phi_center += 2*M_PI;
393 while(phi_center>2*M_PI) phi_center -= 2*M_PI;
394 if(r_center<r_cylinder/2) return -9999;
395 double dphi = acos( r_cylinder/(2*r_center) );
396 if(charge<0) phi = phi_center + dphi;
397 else phi = phi_center - dphi;
398 while(phi<0) phi += 2*M_PI;
399 while(phi>2*M_PI) phi -= 2*M_PI;
400 return phi;
401}
402
403void fill_histogram(double x, double y, double r, int nbin, TH2D* hough)
404{
405 const double M_PI = 3.1415926;
406 double U = 2*x/(x*x+y*y-r*r);
407 double V = 2*y/(x*x+y*y-r*r);
408 double R = 2*r/(x*x+y*y-r*r);
409 double delta_phi = 2.*M_PI/nbin;
410 double phi = -delta_phi/2;
411 for(int i =0; i<nbin; i++){
412 phi += delta_phi;
413 double u = U + R*cos(phi);
414 double v = V + R*sin(phi);
415 double normal = (v-V)/(u-U);
416 double k = -1./normal;
417 double b = v - k*u;
418 double u_cross = -b/(k+1/k);
419 double v_cross = b/(1+k*k);
420
421 double rho=sqrt(u_cross*u_cross+v_cross*v_cross);
422 double theta=atan2(v_cross,u_cross);
423 double slant = V*cos(theta)-U*sin(theta);
424 if( fabs(slant)<0.03 ) continue;
425 if(theta<0) {
426 theta=theta+M_PI;
427 rho=-rho;
428 }
429 if( normal ==0 && x>0) {
430 rho = fabs(x);
431 theta = 0;
432 }
433 if( normal ==0 && x<0) {
434 rho = -fabs(x);
435 theta = M_PI;
436 }
437 if(fabs(rho)>0.1)continue;
438 hough->Fill(theta,rho);
439 }
440}
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
fill_histogram
void fill_histogram(double x, double y, double r, int nbin, TH2D *hough)
Definition: compare_trk_par.cxx:403
intersect_cylinder
double intersect_cylinder(int charge, double x_center, double y_center, double r_cylinder)
Definition: compare_trk_par.cxx:385
compare_trk_par
int compare_trk_par()
Definition: compare_trk_par.cxx:4
t
int t()
Definition: t.c:1