StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
plot.C
1 static const int NQ=4;
2 static const int NL=3;
3 static const int NS=21;
4 static const int NID=NQ*NL*NS;
5 
6 static const int NCUT1=6;
7 static const int NCUT2=6;
8 
9 TH2F* mH2[NQ][NL][NCUT1];
10 TH1F* mHd[NQ][NL][NCUT1];
11 TH2F* mHd2[NQ][NL][NCUT1];
12 
13 static int LOG=0;
14 static TCanvas *c1;
15 static TCanvas *c2;
16 static TCanvas *c3;
17 static TString* FILENAME;
18 TFile* mTFile;
19 
20 //for single point histogram
21 const int col1[NCUT1]={1,2,4,6,8,9};
22 const char* CCUT1[NCUT1]={"All","E>6GeV","gamma","hadron","electron","others"};
23 
24 //for point pair histogram
25 const int col2[NCUT2]={1,2,4,6,8,9};
26 const char* CCUT2[NCUT2]={"All/50","E1,E2>6GeV,E>20GeV","gg","hh","ee","others"};
27 
28 static const float zfms=730.0;
29 static float dx[NID];
30 static float dy[NID];
31 static float dz[NID];
32 static float xx[NID];
33 static float yy[NID];
34 static float zz[NID];
35 
36 static float FIT[NID];
37 static float SIG[NID];
38 static float GEOM[NID];
39 
40 int getSlatid(int q, int l, int s){
41  return (q-1)*NL*NS + (l-1)*NS + s-1;
42 }
43 
44 int getQLS(int slatid, int *q, int *l, int *s){
45  s = slatid%NQ + 1;
46  l = (slatid/NQ)%NL + 1;
47  q = (slatid/NQ/NL)%NQ + 1;
48 }
49 
50 float project(float x, float z, float zp, float vz){
51  //printf("x=%6.2f z=%6.2f zp=%6.2f vtx=%6.2f proj=%6.2f\n",x,z,zp,vz, x/(z-vz)*(zp-vz));
52  return x/(z-vz)*(zp-vz);
53 }
54 
55 void readFpsPosition(char* input="fpsgeom.txt"){
56  FILE *FP = fopen(input,"r");
57  if(!FP) { printf("Could not open %s\n",input); exit;}
58  printf("Reading %s\n",input);
59  char line[1000];
60  int slatid,q,l,s,n=0;
61  float ix,iy,iz,idx,idy,idz;
62  while(fgets(line,1000,FP)!=NULL){
63  sscanf(line,"%d %d %d %d %f %f %f %f %f %f",
64  &slatid,&q,&l,&s,&idx,&idy,&idz,&ix,&iy,&iz);
65  //printf("Id=%3d Q=%3d L=%3d S=%3d D=%9.4f %9.4f %9.4f X=%9.4f %9.4f %9.4f\n",
66  // slatid,q,l,s,idx,idy,idz,ix,iy,iz);
67  xx[slatid]=ix;
68  yy[slatid]=iy;
69  zz[slatid]=iz;
70  dx[slatid]=idx;
71  dy[slatid]=idy;
72  dz[slatid]=idz;
73  //if(q==1 && l==3)
74  if(0)
75  printf("Id=%3d Q=%3d L=%3d S=%3d x=%9.4f %9.4f %9.4f d=%9.4f %9.4f %9.4f\n",
76  slatid,q,l,s,
77  xx[slatid],yy[slatid],zz[slatid],
78  dx[slatid],dy[slatid],dz[slatid]);
79  n++;
80  }
81  printf("Found %d entries\n",n);
82 }
83 
84 void plotd(int quad=0, int layer=0, int cut=0){
85  gStyle->SetOptStat(0);
86  c1->Clear();
87  TVirtualPad *pad;
88  if(quad>0 || layer>0) {pad = c1->cd(0);}
89  for(int l=1; l<=NL; l++){
90  if(layer>0 && layer!=l) continue;
91  if(layer==0){c1->Clear();}
92  if(quad==0) {c1->Divide(2,2);}
93  for(int q=1; q<=NQ; q++){
94  if(quad>0 && quad!=q) continue;
95  if(quad==0) {pad=c1->cd(q);}
96  pad->SetLogy(LOG);
97  mHd[q-1][l-1][cut]->Draw();
98  }
99  if(layer==0){c1->SaveAs(Form("plot2/fmsfpsd_L%1d.png",layer));}
100  }
101  if(quad>0 || layer>0) {c1->SaveAs(Form("plot2/fmsfpsd_Q%1dL%1d.png",quad,layer));}
102 }
103 
104 void plot2d(int quad=0, int layer=0, int slat=0, int cut=0, int val=0, int scale=1, int slice=1){
105  gStyle->SetOptStat(0);
106  c1->Clear();
107  TVirtualPad *pad;
108  if(quad>0 || layer>0) {pad = c1->cd(0);}
109  for(int l=1; l<=NL; l++){
110  if(layer>0 && layer!=l) continue;
111  if(layer==0){c1->Clear();}
112  if(slice==0 && quad==0) {c1->Divide(2,2);}
113  for(int q=1; q<=NQ; q++){
114  if(quad>0 && quad!=q) continue;
115  if(quad==0) {pad=c1->cd(q);}
116  pad->SetLogz(LOG);
117  TH2F* h;
118  if(val==0) h=mH2[q-1][l-1][cut];
119  if(val==1) h=mHd2[q-1][l-1][cut];
120  TH1D* projx=h->ProjectionX();
121  TH2F* h2=(TH2F*)h->Clone(Form("Q%1dL%1d_scaled",q,l)); h2->Reset();
122  int nx=h->GetNbinsX();
123  int ny=h->GetNbinsY();
124  for(int x=0; x<=nx+1; x++){
125  float w=projx->GetBinContent(x);
126  for(int y=0; y<=ny+1; y++){
127  int bin=h->GetBin(x,y);
128  float v=h->GetBinContent(bin);
129  float newv=0.0;
130  if(w>0.0) newv=v/w;
131  h2->SetBinContent(bin,newv);
132  //printf("x=%3d y=%3d w=%6.2f v=%6.2f v/w=%8.5f\n",x,y,w,v,newv);
133  }
134  }
135  if(slice==0){
136  //h2->SetMinimum(h2->GetMaximum()*0.6);
137  if(scale) {h2->Draw("colz");}
138  else {h->Draw("colz");}
139  }else{
140  c1->Clear();
141  if(slat==0) {c1->Divide(3,7);}
142  for(int s=1; s<=NS; s++){
143  int slatid=getSlatid(q,l,s);
144  if(slat>0 && s!=slat) continue;
145  if(s>19 && (q==2 || q==4)) continue;
146  char cc[100];
147  sprintf(cc,"Q%1dL%1dS%02d",q,l,s);
148  TH1D *hp = h2->ProjectionX(cc,s,s);
149  float xxx,dxx;
150 
151  if(0){ //playing by hand
152  if(l==1){
153  xxx=xx[slatid]; dxx=dx[slatid];
154  if(q<=2) xxx=xxx/1.02 - 0.5;
155  if(q>=3) xxx=xxx*1.015 + 1.5;
156  }else if(l==2){
157  xxx=yy[slatid]; dxx=dy[slatid];
158  xxx=xxx/1.010 + 0.5;
159  }else if(l==3){
160  xxx=yy[slatid]; dxx=dy[slatid];
161  xxx=xxx/1.018 + 0.5;
162  }
163  }
164  if(1){ //offsets are now in geom file
165  float factor=1.00;
166  if(l==1){
167  xxx=xx[slatid]; dxx=dx[slatid];
168  if(q<=2) xxx=xxx/factor;
169  //if(q>=3) xxx=xxx*1.015;
170  if(q>=3) xxx=xxx/factor;
171  }else if(l==2){
172  xxx=yy[slatid]; dxx=dy[slatid];
173  xxx=xxx/factor;
174  }else if(l==3){
175  xxx=yy[slatid]; dxx=dy[slatid];
176  xxx=xxx/factor;
177  }
178  }
179 
180  float xmin=xxx-3.0*dxx;
181  float xmax=xxx+3.0*dxx;
182  printf("xmin/xmax=%f %f\n",xmin,xmax);
183  int bmin=hp->FindBin(xmin);
184  int bmax=hp->FindBin(xmax);
185  printf("bmin/bmax=%f %f\n",bmin,bmax);
186  hp->GetXaxis()->SetRange(bmin,bmax);
187  float min=hp->GetMinimum(0.02);
188  float max=hp->GetMaximum();
189  float delta=max-min;
190  float peak=hp->GetBinCenter(hp->GetMaximumBin());
191  float avg=hp->Integral(bmin,bmax)/(bmax-bmin+1);
192  float sig=dxx/2.0/2.0;
193  //printf("min/max/delta/peak/avg/sig=%f %f %f %f %f %f\n",min,max,delta,peak,avg,sig);
194  float ymin=min-delta*0.1;
195  float ymax=max+delta*0.3;
196  //hp->SetMaximum(ymax);
197  //hp->SetMinimum(ymin);
198  if(slat==0){
199  pad=c1->cd(s);
200  float mergin=0.005;
201  pad->SetRightMargin(mergin); pad->SetLeftMargin(mergin);
202  pad->SetTopMargin(mergin); pad->SetBottomMargin(0.01);
203  }
204  if(q==1 || quad>0 || slat>0) {hp->Draw();}
205  else {hp->Draw("same");}
206  //draw lines from geometry parameters
207  for(int i=-1; i<2; i++){
208  float proj=xxx+i*dxx*0.5;
209  if(i==0) GEOM[slatid]=proj;
210  //printf("proj,ymin,ymax=%f %f %f\n",proj,ymin,ymax);
211  TLine *ll = new TLine(proj,ymin,proj,ymax);
212  ll->SetLineColor(kBlue); ll->SetLineWidth(2);
213  if(i!=0) ll->Draw();
214  }
215  //fit
216  TF1 *f=new TF1("ff", "gaus(0)+[3]", xmin, xmax);
217  xmin=proj-3.0*dxx;
218  xmax=proj+3.0*dxx;
219  //printf("xmin=%6.3f xmax=%6.3f max=%6.3f avg=%6.3f max-avg=%6.3f geom=%6.3f dxx=%6.3f\n",
220  // xmin,xmax,max,avg,max-avg,GEOM[slatid],dxx);
221  f->SetParameter(0,max-avg);
222  //f->SetParameter(1,GEOM[slatid]);
223  f->SetParameter(1,peak);
224  f->SetParameter(2,sig);
225  f->SetParameter(3,avg);
226  int res = hp->Fit("ff","0R");
227  FIT[slatid]=f->GetParameter(1);
228  SIG[slatid]=f->GetParameter(2);
229  //printf("A=%6.3f FIT=%6.3f SIG=%6.2f Base=%6.3f\n",f->GetParameter(0),FIT[slatid],SIG[slatid],f->GetParameter(3));
230  f->SetLineColor(kRed);
231  f->Draw("same");
232  } //loop over slat
233  c1->SaveAs(Form("plot2/fmsfps2_Q%1dL%1d_cut%d_slice.png",q,layer,cut));
234  } //if slice==1
235  } //loop over quad
236  if(slice==0) c1->SaveAs(Form("plot2/fmsfps2_Q%1dL%1d_cut%d.png",quad,layer,cut));
237  } //loop over layer
238  //if(quad>0 || layer>0) {c1->SaveAs(Form("plot2/fmsfps2_Q%1dL%1d_cut%d.png",quad,layer,cut));}
239  if(slice==1 && slat==0){
240  for(int l=1; l<=NL; l++){
241  if(layer>0 && l!=layer) continue;
242  c2->Clear();
243  if(quad==0) c2->Divide(2,2);
244  for(int q=1; q<=NQ; q++){
245  if(quad>0 && q!=quad) continue;
246  if(quad==0) c2->cd(q);
247  TGraphErrors *g=new TGraphErrors(1);
248  char t[100];
249  sprintf(t,"Q%1dL%1d Fit-DB vs slat",q,l);
250  g->SetTitle(t);
251  float x1,x2;
252  int pm[NQ][NL]={{1,1,1},{1,-1,-1},{-1,1,1},{-1,-1,-1}};
253  if(pm[q-1][l-1]==1) {x1=-10.0; x2=110.0;}
254  else {x1=-110.0; x2=10.0;}
255  TH2F *frame=new TH2F(t,t,1,x1,x2,1,-10.0,10.0);
256  frame->Draw();
257  int np=0;
258  for(int s=1; s<=NS; s++){
259  int slatid=getSlatid(q,l,s);
260  double d=FIT[slatid]-GEOM[slatid];
261  if(fabs(d)<2.0 && SIG[slatid]>0.5 && SIG[slatid]<6.0){
262  g->SetPoint(np,GEOM[slatid],d);
263  g->SetPointError(np,0.0,SIG[slatid]);
264  np++;
265  }
266  }
267  g->SetMarkerStyle(20); g->SetMarkerColor(kRed);
268  g->Draw("pw");
269  //fit
270  TF1 *f2=new TF1("f2", "[0]+[1]*x", -100, 100);
271  f2->SetParameter(0,0.0);
272  f2->SetParameter(1,0.0);
273  int res = g->Fit("f2","0RQ");
274  f2->SetLineColor(kBlue);
275  f2->Draw("same");
276  }
277  c2->SaveAs(Form("plot2/fmsfps2_Q%1dL%1d_cut%d_align.png",quad,layer,cut));
278  }
279  }
280 }
281 
282 void pltcut(char* name, TCanvas* c, int div, int log=0, int legend=0){
283  char cc[100];
284  TString hname(name);
285  if(log==0) c->cd(div);
286  if(log==1) c->cd(div)->SetLogy();
287  if(name!=""){
288  for(int cut=0; cut<NCUT1; cut++){
289  sprintf(cc,"%s_c%d",name,cut);
290  TH1* h=(TH1*)mTFile->Get(cc);
291  if(cut==0 && hname.Contains("p_")) h->Scale(0.05);
292  h->SetLineColor(col1[cut]); h->SetMarkerColor(col1[cut]);
293  if(log==1) h->SetMinimum(0.1);
294  if(name=="NP") h->SetMaximum(h->GetMaximum()*10.0);
295  TString opt = "";
296  if(cut>0) opt+="same";
297  h->Draw(opt.Data());
298  }
299  }
300  if(legend){
301  for(int cut=0; cut<NCUT1; cut++){
302  TText* t=new TText(0.7,0.6-0.07*cut,CCUT1[cut]);
303  t->SetTextSize(0.07); t->SetNDC(); t->SetTextColor(col1[cut]);
304  t->Draw();
305  }
306  }
307 }
308 
309 void pltcut2(char* name, TCanvas* c, int div, int log=0, int legend=0){
310  char cc[100];
311  TString hname(name);
312  if(log==0) c->cd(div);
313  if(log==1) c->cd(div)->SetLogy();
314  if(log==2) c->cd(div)->SetLogz();
315  if(name!=""){
316  for(int cut=0; cut<NCUT2; cut++){
317  if(name=="p_pid" && cut==0) continue;
318  sprintf(cc,"%s_c%d",name,cut);
319  TH1* h=(TH1*)mTFile->Get(cc);
320  if(cut==0) h->Scale(1/40.0);
321  if(cut==3 && name=="p_m1") h->Scale(20.0);
322  if(cut==4 && name=="p_m1") h->Scale(40.0);
323  h->SetLineColor(col2[cut]); h->SetMarkerColor(col2[cut]);
324  if(log==1) h->SetMinimum(0.1);
325  if(name=="p_NP") h->SetMaximum(h->GetMaximum()*1000.0);
326  TString opt = "";
327  if(cut>0) opt+="same";
328  //if(name=="p_ept") opt+="box";
329  if(name=="p_pid") opt+="colz";
330  h->Draw(opt.Data());
331  }
332  }
333  if(legend){
334  for(int cut=0; cut<NCUT2; cut++){
335  TText* t=new TText(0.7,0.6-0.07*cut,CCUT2[cut]);
336  t->SetTextSize(0.07); t->SetNDC(); t->SetTextColor(col2[cut]);
337  t->Draw();
338  }
339  }
340 }
341 
342 void plotfms(){
343  c1->Clear();
344  c1->Divide(2,3);
345  pltcut("NP", c1,1,1,1);
346  pltcut("e", c1,2,1);
347  pltcut("elo",c1,3,1);
348  pltcut("pt", c1,4,1);
349  pltcut("ept",c1,5,0);
350  pltcut("eta",c1,6,1);
351  c1->SaveAs("plot2/fms1.png");
352 
353  c2->Clear();
354  c2->Divide(2,3);
355  pltcut("phi", c2,1,1,1);
356  pltcut("x", c2,2,1);
357  pltcut("y", c2,3,1);
358  pltcut("xy", c2,4,0);
359  pltcut("pid", c2,5,0);
360  pltcut("pid2",c2,6,0);
361  c2->SaveAs("plot2/fms2.png");
362 
363  c3->Clear();
364  c3->Divide(2,3);
365  c3->cd(1); xy_c0->Draw("");
366  c3->cd(2); xy_c1->Draw("");
367  c3->cd(3); xy_c2->Draw("");
368  c3->cd(4); xy_c3->Draw("");
369  c3->cd(5); xy_c4->Draw("");
370  c3->cd(6); xy_c5->Draw("");
371  c2->SaveAs("plot2/fms3.png");
372 }
373 
374 void plotfmsPair(){
375  c1->Clear();
376  c1->Divide(2,3);
377  pltcut2("p_NP", c1, 1, 1);
378  pltcut2("p_e", c1, 2, 1);
379  pltcut2("p_pt", c1, 3, 1);
380  pltcut2("p_ept", c1, 4, 0);
381  pltcut2("p_eta", c1, 5, 1);
382  pltcut2("p_phi", c1, 6, 1, 1);
383  c1->SaveAs("plot2/pair1.png");
384 
385  c2->Clear();
386  c2->Divide(2,3);
387  pltcut2("p_m1", c2, 1, 0);
388  pltcut2("p_m2", c2, 2, 1);
389  pltcut2("p_zgg", c2, 3, 1);
390  pltcut2("p_r30", c2, 4, 1);
391  pltcut2("p_r100", c2, 5, 1, 1);
392  pltcut2("p_pid", c2, 6, 2);
393  c2->SaveAs("plotpair2.png");
394 }
395 
396 void readfile(int file=0){
397  switch(file){
398  case 0: mTFile = new TFile("hist/fmsps.root","old"); break;
399  case 1: mTFile = new TFile("sim/test_electron.fmsfps.root"); break;
400  case 2: mTFile = new TFile("sim/test_pion.fmsfps.root"); break;
401  case 3: mTFile = new TFile("sim/test_gamma.fmsfps.root"); break;
402  case 4: mTFile = new TFile("sim/test_pi0.fmsfps.root"); break;
403  case 5: mTFile = new TFile("sim/test_mu-.fmsfps.root"); break;
404  }
405  char c[100];
406  for(int cut=0; cut<NCUT1; cut++){
407  for(int q=0; q<NQ; q++){
408  for(int l=0; l<NL; l++){
409  sprintf(c,"FMSFPS_Q%1dL%1d_c%d",q+1,l+1,cut);
410  mH2[q][l][cut]=(TH2F*)mTFile->Get(c);
411  sprintf(c,"FMS-FPS_Q%1dL%1d_c%d",q+1,l+1,cut);
412  mHd[q][l][cut]=(TH1F*)mTFile->Get(c);
413  sprintf(c,"FMSFPSd_Q%1dL%1d_c%d",q+1,l+1,cut);
414  mHd2[q][l][cut]=(TH2F*)mTFile->Get(c);
415  }
416  }
417  }
418 }
419 
420 void openCanvas(){
421  c1 = new TCanvas("FPS","FPS",50,10,700,800);
422  c2 = new TCanvas("FPS2","FPS2",750,10,700,800);
423  c3 = new TCanvas("FPS3","FPS3",1450,10,700,800);
424  gStyle->SetPalette(1);
425  gStyle->SetStatW(0.4);
426 }
427 
428 void plot(int file=0, int plt=0, int quad=0, int layer=0, int slat=0, int cut=0, int log=0){
429  LOG=log;
430  readFpsPosition();
431  readfile(file);
432  openCanvas();
433  if(plt==1 || plt==0) plotfms();
434  if(plt==2 || plt==0) plotfmsPair();
435  if(plt==3 || plt==0) plot2d(quad,layer,slat,cut,0,1,0);
436  if(plt==4 || plt==0) plot2d(quad,layer,slat,cut,0,1,1);
437  if(plt==5 || plt==0){
438  for(int l=1; l<=NL; l++){
439  plot2d(0,l,0,cut,0,1,0);
440  plot2d(0,l,0,cut,0,1,1);
441  }
442  }
443 }
Generated by   doxygen 1.8.5