StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
UtilBeamLine3D.cxx
1 #include <math.h>
2 #include <assert.h>
3 #include <TObjArray.h>
4 #include <TH1.h>
5 #include <TH2.h>
6 #include <TLine.h>
7 #include <TList.h>
8 
9 #include "UtilBeamLine3D.h"
10 
11 //========================
12 //========================
13 UtilBeamLine3D::UtilBeamLine3D() {
14  track.clear(); // just on case;
15  fcnCount=0;
16  // input tracks
17  cut_maxRxy=1.5;
18  cut_maxZ=100;
19  cut_minChi2=0.4;
20  cut_maxChi2=2.;
21  cut_minP=1.2; // GeV/c
22  // minimization:
23  cut_Dmax=1.0; // likelhood cut-off
24  par_filter=0; // 0=all, 1: West-only = +Z & +eta
25  printf("Params of UtilBeamLine3D\n INPUT tracks: maxRxy=%.1fcm, maxZ=%.1fcm, Chi2=[%.1f,%.1f] minP=%.1f GeVC/c \n Minimization: Dmax=%.1fcm, WestOnly=%d \n\n",
26  cut_maxRxy,cut_maxZ,cut_minChi2, cut_maxChi2,cut_minP,cut_Dmax,par_filter);
27 }
28 
29 //========================
30 //========================
31 void
32 UtilBeamLine3D:: initHisto( TObjArray * HList){
33  TList *L=0; TLine *ln=0; float yMax=1e6; TH1* h=0;
34 
35  memset(hA,0,sizeof(hA));
36  printf("helper:initHist\n");
37 
38  hA[0]=h=new TH1F("trStat","statistics for tracks",20,0.5,20.5);
39  h->GetXaxis()->SetTitleOffset(0.4); h->GetXaxis()->SetLabelSize(0.06); h->GetXaxis()->SetTitleSize(0.05); h->SetMinimum(0.8);
40  h->SetLineColor(kBlue);h->SetLineWidth(2);
41  h->SetMarkerSize(2);//<-- large text
42 
43  hA[1]=h=new TH1F("trCh2","input track chi2; chi2/dof",100,0,10);
44  ln=new TLine(cut_minChi2,0,cut_minChi2,yMax);
45  ln->SetLineColor(kRed); ln->SetLineStyle(2);
46  L= h->GetListOfFunctions(); L->Add(ln);
47  ln=new TLine(cut_maxChi2,0,cut_maxChi2,yMax);
48  ln->SetLineColor(kRed); ln->SetLineStyle(2);
49  L= h->GetListOfFunctions(); L->Add(ln);
50 
51  hA[2]=h=new TH1F("trZ","input track Z @DCA; Z(cm)",50,-200,200);
52  ln=new TLine(cut_maxZ,0,cut_maxZ,yMax);
53  ln->SetLineColor(kRed); ln->SetLineStyle(2);
54  L= h->GetListOfFunctions(); L->Add(ln);
55  ln=new TLine(-cut_maxZ,0,-cut_maxZ,yMax);
56  ln->SetLineColor(kRed); ln->SetLineStyle(2);
57  L= h->GetListOfFunctions(); L->Add(ln);
58 
59  hA[3]=h=new TH1F("trX","input track X @DCA; X(cm)",100,-4,4);
60  ln=new TLine(cut_maxRxy,0,cut_maxRxy,yMax);
61  ln->SetLineColor(kRed); ln->SetLineStyle(2);
62  L= h->GetListOfFunctions(); L->Add(ln);
63  ln=new TLine(-cut_maxRxy,0,-cut_maxRxy,yMax);
64  ln->SetLineColor(kRed); ln->SetLineStyle(2);
65  L= h->GetListOfFunctions(); L->Add(ln);
66 
67  hA[4]=h=new TH1F("trY","input track Y @DCA; Y(cm)",100,-4,4);
68  ln=new TLine(cut_maxRxy,0,cut_maxRxy,yMax);
69  ln->SetLineColor(kRed); ln->SetLineStyle(2);
70  L= h->GetListOfFunctions(); L->Add(ln);
71  ln=new TLine(-cut_maxRxy,0,-cut_maxRxy,yMax);
72  ln->SetLineColor(kRed); ln->SetLineStyle(2);
73  L= h->GetListOfFunctions(); L->Add(ln);
74 
75  hA[5]=h=new TH1F("trP","input track momentum P @DCA; P (GeV/c)",100,0,10);
76  ln=new TLine(cut_minP,0,cut_minP,yMax);
77  ln->SetLineColor(kRed); ln->SetLineStyle(2);
78  L= h->GetListOfFunctions(); L->Add(ln);
79 
80  hA[6]=new TH1F("trPt","input track momentum PT @DCA; PT (GeV/c)",100,0,10);
81  //free 7..12
82 
83  int nB=30;
84  hA[13]=h= new TH2D("fcnChiXY","ln(3D likelihood), STAR X-Y plane; X-axis (cm); Y-axis(cm)",nB,-.3,.3,nB,-.3,.3);
85  ln=new TLine(-10,0,10,0); ln->SetLineColor(kMagenta);
86  L= h->GetListOfFunctions(); L->Add(ln);
87  ln=new TLine(0,-10,0,10); ln->SetLineColor(kMagenta);
88  L= h->GetListOfFunctions(); L->Add(ln);
89 
90  hA[14]=h= new TH2D("fcnChiXnX","ln(3D likelihood), STAR X-n_{X} plane; X-axis (cm); slope X ",nB,0,1,nB,-2e-3,2e-3);
91  ln=new TLine(0,-10,0,10); ln->SetLineColor(kMagenta);
92  L= h->GetListOfFunctions(); L->Add(ln);
93  ln=new TLine(-10,0,10,0); ln->SetLineColor(kMagenta);
94  L= h->GetListOfFunctions(); L->Add(ln);
95 
96  hA[15]=h= new TH2D("fcnChiYnY","ln(3D likelihood), STAR Y-n_{Y} plane; Y-axis (cm); slope Y ",nB,0,1,nB,-2e-3,2e-3);
97  ln=new TLine(0,-10,0,10); ln->SetLineColor(kMagenta);
98  L= h->GetListOfFunctions(); L->Add(ln);
99  ln=new TLine(-10,0,10,0); ln->SetLineColor(kMagenta);
100  L= h->GetListOfFunctions(); L->Add(ln);
101  // free 16..18
102 
103  hA[19]=h=new TH1F("bmSol","Beam line params (3D fit); value",4,1,1);
104  h->SetMarkerStyle(24);
105  float mY=0.3;
106  h->SetMaximum(mY); h->SetMinimum(-mY);
107  h->GetXaxis()->SetTitleOffset(0.4); h->GetXaxis()->SetLabelSize(0.06); h->GetXaxis()->SetTitleSize(0.05);
108  h->SetLineColor(kBlue);h->SetLineWidth(2);
109  // h->SetMarkerSize(2);//<-- large text
110 
111  ln=new TLine(0,0,6,0); L= h->GetListOfFunctions(); L->Add(ln);
112 
113  //20-30: reserved for chi2 monitoring
114  hA[20]=new TH1F("fcnDet","FCN: determinant from DCA solution",200,0,1.01);
115  hA[21]=h=new TH1F("fcnDca","FCN: 3D dca length ...; DCA (cm)",100,0,5.);
116  ln=new TLine(cut_Dmax,0,cut_Dmax,yMax);
117  ln->SetLineColor(kRed); ln->SetLineStyle(2);
118  L= h->GetListOfFunctions(); L->Add(ln);
119 
120  hA[22]=h=new TH2F("fcnDcaZ","FCN: 3D dca length; Z (cm);DCA (cm)",50,-200,200.,50,0,5.);
121  ln=new TLine(-200,cut_Dmax,200,cut_Dmax);
122  ln->SetLineColor(kRed); ln->SetLineStyle(2);
123  L= h->GetListOfFunctions(); L->Add(ln);
124  hA[23]=h=new TH1F("fcnNtr","FCN: # of tracks with 3D.DCA<Dmax ; Z (cm)",50,-200,200.);
125  h->SetLineColor(kBlue); h->SetFillColor(kBlue);
126 
127  for(int i=0;i<mxH; i++) if(hA[i]) HList->Add(hA[i]);
128  printf("track QA params: maxRxy<%.1f cm, maxZ<%.1f cm, chi2/dof=[%.1f, %.1f], Dmax=%.2f cm track filter=%d\n", cut_maxRxy, cut_maxZ, cut_minChi2,cut_maxChi2,cut_Dmax,par_filter);
129 
130 }
131 
132 
133 //========================
134 //========================
135 int
136 UtilBeamLine3D::qaTracks(){
137  printf("qaTRacks start=%d\n",(int)track.size());
138  int nOK=0;
139  vector<TrackStump>::iterator it;
140  for (it=track.begin() ; it < track.end(); it++ ) {
141  TrackStump *t= &(*it);
142  hA[0]->Fill(1);
143 
144  hA[1]->Fill(t->chi2);
145  if( t->chi2<cut_minChi2 || t->chi2>cut_maxChi2) continue;
146  hA[0]->Fill(2);
147 
148  hA[2]->Fill(t->r.z());
149  if(fabs(t->r.z()) >cut_maxZ) continue;
150  hA[0]->Fill(3);
151 
152  hA[3]->Fill(t->r.x());
153  if(fabs(t->r.x()) >cut_maxRxy) continue;
154  hA[0]->Fill(4);
155 
156  hA[4]->Fill(t->r.y());
157  if(fabs(t->r.y()) >cut_maxRxy) continue;
158  hA[0]->Fill(5);
159 
160  hA[5]->Fill(t->P);
161  if(t->P<cut_minP) continue;
162  hA[0]->Fill(6);
163 
164  if(par_filter==1) {
165  if(t->r.z()< 0 || t->p.z() <0) continue;
166  }
167  hA[0]->Fill(7);
168 
169 
170  //..... track has been accepted ....
171  t->bad=0;
172  nOK++;
173  hA[6]->Fill(t->Pt);
174  }
175  printf("qaTRacks end=%d\n",nOK);
176  return nOK;
177 }
178 
179 
180 //========================
181 //========================
182 void
183 UtilBeamLine3D::print(int k) {
184  printf("UtilBeamLine3D::print(%d) size=%d\n",k,(int)track.size());
185  for(int i=0;i<(int)track.size();i++) {
186  track[i].print();
187  if(i>=k) break;
188  }
189 }
190 
191 //========================
192 //========================
193 void
194 UtilBeamLine3D::scanChi2(double *par, int mode){
195 
196  // mode: //0=x-y, 1=x-nx , 2= y-ny
197  printf("scan chi2, mode=%d ....\n",mode);
198  // working variables for likelihod fcn
199  int npar,iflag=0;
200  double *grad=0;
201  double fcnval;
202  double wrkPar[mxPar];
203 
204  // 2D scan in X-Y plane
205  // copy fixed Ux,Uy
206 
207  TH2F *h2=0;
208  if(mode==0) h2=(TH2F*)hA[13];
209  if(mode==1) h2=(TH2F*)hA[14];
210  if(mode==2) h2=(TH2F*)hA[15];
211  TAxis *xax=h2->GetXaxis();
212  TAxis *yax=h2->GetYaxis();
213 
214  if(mode==0) { // adjust only if too much displaced
215  if(fabs(par[0])>0.2 ||fabs(par[1])>0.2 ){
216  xax->Set(xax->GetNbins(),par[0]-0.3,par[0]+0.3);
217  yax->Set(yax->GetNbins(),par[1]-0.3,par[1]+0.3);
218  }
219  }
220  if(mode>0) { // always adjust
221  xax->Set(xax->GetNbins(),par[0]-0.3,par[0]+0.3);
222  yax->Set(yax->GetNbins(),par[2]-9e-3,par[2]+9e-3);
223  }
224 
225  // fill chi2 values
226  for(int bx=1;bx<=xax->GetNbins();bx++){
227  double x=xax->GetBinCenter(bx);
228  for(int by=1;by<=yax->GetNbins();by++){
229  double y=yax->GetBinCenter(by);
230  if(mode==0) { // X-Y plane
231  wrkPar[0]=x; wrkPar[1]=y;
232  wrkPar[2]=par[2]; wrkPar[3]=par[3];
233  }
234  if(mode==1) { // X-nX
235  wrkPar[0]=x; wrkPar[1]=par[1];
236  wrkPar[2]=y; wrkPar[3]=par[3];
237  }
238 
239  if(mode==2) { // X-nX
240  wrkPar[0]=par[0]; wrkPar[1]=x;
241  wrkPar[2]=par[2]; wrkPar[3]=y;
242  }
243 
244  beamLineLike3D(npar,grad,fcnval, wrkPar,iflag);
245  //printf(" %d %f %f fcn=%.1f\n",by,x,y,fcnval);
246  h2->Fill(x,y,fcnval);
247  }
248  }
249 
250 }
251 
252 //========================
253 void
254 UtilBeamLine3D::evalSolution(double *par){
255  // working variables for likelihod fcn
256  int npar,iflag=0;
257  double *grad=0;
258  double fcnval;
259  // double wrkPar[mxPar];
260  fcnMon1=1; // activate chi2 monitor for the 1st pass
261  beamLineLike3D(npar,grad,fcnval, par,iflag);
262 }
263 
264 
265 //========================
266 //========================
267 void
268 UtilBeamLine3D::readTracks(const TString fnameT){
269  const char *fname=fnameT.Data();
270  printf("Read fname=%s=\n",fname);
271  FILE *fd=fopen(fname,"r"); assert(fd);
272  char text[1000];
273  TrackStump t;
274  float x,y,z,px,py,pz,ery,eryz,erz;
275  track.clear();
276  int mxSkip=2; // max # consecutive of lines allowed to be wrong
277  int nSkip=mxSkip;
278  char buf[3000];
279  while(1) {
280  char * retc=fgets(buf,3000,fd);
281  if(retc==0) break;
282  //if(buf[0]=='#') continue;
283  //printf("=%s=\n",buf);
284  int ret=sscanf(buf ,"%s %f %f %f %f %f %f %f %f %f %d %f %f %d,",text,&x,&y,&z,&px,&py,&pz,&ery,&eryz,&erz,&t.nFitP,&t.chi2,&t.z0,&t.eveId);
285  //printf("%d %s %f %f %f %f %f %f %f %f %f %d %f %.1f %d \n",ret,text,x,y,z,px,py,pz,t.cyy,t.czy,t.czz,t.nFitP,t.chi2,t.z0,t.eveId);
286  if(ret!=14) {
287  printf("bad ret=%d nSkip=%d size=%d\n",ret,nSkip,(int)track.size());
288  nSkip--;
289  if(nSkip) continue; // allow this line as bad ;
290  else break; // do not read from input any more
291  }
292  nSkip=mxSkip;
293  // add gentle slope in x & y
294  //x=0.0002*z; y=-0.0003*z;
295  t.r.SetXYZ(x,y,z);
296  t.p.SetXYZ(px,py,pz);
297  t.P=t.p.Mag();
298  t.Pt=t.p.Pt();
299  t.p=t.p.Unit();
300  t.ery2=ery*ery;
301  t.eryz=eryz;
302  t.erz2=erz*erz;
303  t.bad=1; // assume all bad, track QA will decide
304  track.push_back(t);
305  //if(track.size()>19) break;
306  }
307  printf("found %d tracks\n",(int)track.size());
308  fclose(fd);
309 }
Generated by   doxygen 1.8.5