StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Vertex3D.cxx
1 #include <string.h>
2 #include <assert.h>
3 #include <StMessMgr.h>
4 
5 #include <TH1.h>
6 #include <TH2.h>
7 #include <TObjArray.h>
8 #include <TList.h>
9 #include <TLine.h>
10 #include <TArrow.h>
11 #include <TEllipse.h>
12 
13 #include <StEvent/StGlobalTrack.h>
14 
15 #include <StBFChain/StBFChain.h> // for geant vertex
16 #include <tables/St_g2t_vertex_Table.h> //for geant vertex
17 
18 #include "Vertex3D.h"
19 
20 namespace StEvPPV {
21 
22 //==========================================================
23 //==========================================================
24 Vertex3D::Vertex3D() {
25 }
26 
27 //==========================================================
28 //==========================================================
29 void
30 Vertex3D::initRun(){
31  LOG_INFO <<Form("Vertex3D::initRun() params:track glob %.1f<PT(GeV/c)<%.1f, sigY<%.1fcm, nTrack>%d ", cut_pT1, cut_pT2, cut_sigY, cut_numTrack)<<endm;
32  }
33 
34 //==========================================================
35 //==========================================================
36 Vertex3D::~Vertex3D(){
37 }
38 
39 //==========================================================
40 //==========================================================
41 void
42 Vertex3D::clearEvent(){
43  hA[0]->Fill(1);
44  isFound=false;
45  track.clear();
46  // printf("ver3DclearEvt\n");
47 }
48 
49 
50 //==========================================================
51 //==========================================================
52 void
53 Vertex3D::clearTracks(){
54  hA[0]->Fill(2);
55  track.clear();
56  // printf("ver3DclearTrk\n");
57 }
58 
59 
60 //==========================================================
61 //==========================================================
62 void Vertex3D::addTrack(TrackData* trk)
63 {
64  if(isValid()) return; // block adding new track after 1st vertex is found
65  if( !trk->mTpc && !trk->mBtof && !trk->mCtb && !trk->mBemc && trk->mEemc ) return;
66  hA[0]->Fill(3);
67  // track is matched to any fast detector
68 
69  float pt=trk->dcaTrack.gP.Pt();
70  hA[2]->Fill(pt);
71  if (pt <cut_pT1) return;
72  hA[0]->Fill(4);
73  if (pt >cut_pT2) return;
74  hA[0]->Fill(5);
75 
76  float sY=trk->dcaTrack.sigYloc; //track extrapolation error in transverse direction
77  hA[3]->Fill(sY);
78  if(sY>cut_sigY) return;
79  hA[0]->Fill(6);
80 
81  hA[4]->Fill(trk->dcaTrack.sigZ);
82  track.push_back(trk);
83 }
84 
85 //==========================================================
86 //==========================================================
87 void Vertex3D::doExtrapolation()
88 { // study track cov matrix for individual tracks using my macro plTrCov.C
89 #if 0 // Looks like this function is doing nothing
90  for(unsigned int i=0;i<track.size();i++) {
91  const StiKalmanTrack* tr=track[i]->mother;
92  hA[5]->Fill(1);
93  if( fabs(tr->getChi2()-1.)>0.8) continue;
94  hA[5]->Fill(2);
95  if( tr->getFitPointCount()<20) continue;
96  hA[5]->Fill(3);
97  if( fabs(tr->getPseudoRapidity())>0.5) continue;
98  hA[5]->Fill(4);
99  if( fabs(tr->getPt())<1.) continue;
100  hA[5]->Fill(5);
101 
102  printf("TTT global phi/rad=%.3f PT=%.2f curv=%f Rxy=%.1f nFitP=%d eta=%.2f chi2/dof=%.1f\n",tr->getPhi(),tr->getPt(),tr->getCurvature(),tr->getPoint().perp(),tr->getFitPointCount(),tr->getPseudoRapidity(),tr->getChi2());
103  printf("T: track position along nodes, N=%d\n",tr->getFitPointCount());
104  StiKTNIterator tNode = tr->rbegin();
105  StiKTNIterator eNode = tr->rend();
106  int iN=0;
107  float minLx=9999, maxLx=0.; // local coordinates of track
108  for (;tNode!=eNode;++tNode){
109  StiKalmanTrackNode *node = &(*tNode);
110  if(!node->isValid()) continue;
111 
112  StiHit *stiHit = node->getHit();
113  if (!stiHit) continue;
114 
115  if (node->getChi2()>1000) continue;
116  if (!node->isFitted()) continue;
117  iN++;
118  float Lx=node->x();
119  float Ly=node->y();
120  float Lz=node->z();
121  float sLy=sqrt(node->fitErrs()._cYY);
122  float sLz=sqrt(node->fitErrs()._cZZ);
123  float Rxy=sqrt(Lx*Lx+Ly*Ly);
124  if(minLx>Lx) minLx=Lx;
125  if(maxLx<Lx) maxLx=Lx;
126  printf("iN=%d Lx=%.2f Ly=%.2f %.2f Lz=%.2f %.2f Rxy=%.1f gX=%.2f gY=%.2f \n",iN,Lx,Ly,sLy,Lz,sLz,Rxy,node->x_g(),node->y_g());
127  printf("TTm %.2f %.2f %.3f %.2f %.3f %f\n",node->x_g(),node->y_g(),sLy,node->z_g(),sLz,node->getAlpha());
128 
129  } // end of nodes w/ hits
130 
131  for(int iDir=0;iDir<2;iDir++) {
132  if(iDir==0) printf("extrapolate track toward vertex, innLx=%.1f\n",minLx);
133  if(iDir==1) printf("extrapolate track toward BSMD, outLx=%.1f\n",maxLx);
134  for(int iI=0;iI<17;iI++) {
135  float Lx=100.;StiKalmanTrackNode* node=0; int ret=-2;
136  StiKalmanTrack track2=*tr; // clone track to allow extrapolation
137  if(iDir==0) {// toward vertex
138  Lx=minLx-iI*5;
139  node=track2.getInnerMostNode();
140  ret=node->propagate(Lx,kCylindrical,kOutsideIn);
141  } else {
142  Lx=maxLx+iI*2; if(Lx>250) break;
143  node=track2.getOuterMostNode();
144  ret=node->propagate(Lx,kCylindrical,kInsideOut);
145  }
146  if(iI==0)cout<<"TT track extrap node:"<< *node;
147  if(ret<0) {printf("prop in x=%f ret=%d, skip\n",Lx,ret); continue;}
148  float Ly=node->y();
149  float Lz=node->z();
150  if(fabs(Ly)>40.) break; // to not cross sector boundary ??
151  if(iDir==1 && fabs(Ly)>36.) break; // to not cross sector boundary ??
152  node->propagateError();
153  float sLy=sqrt(node->fitErrs()._cYY);
154  float sLz=sqrt(node->fitErrs()._cZZ);
155  printf("%diN=%d Lx=%.2f Ly=%.2f %.2f Lz=%.2f %.2f gX=%.2f gY=%.2f\n",iDir,iI,Lx,Ly,sLy,Lz,sLz,node->x_g(),node->y_g());
156  printf("TT%d %.2f %.2f %.3f %.2f %.3f %f\n",iDir,node->x_g(),node->y_g(),sLy,node->z_g(),sLz,node->getAlpha());
157  }
158  }
159 
160  }
161 #endif
162 }
163 //head -n 4224 Lda83-618 | tail -n 150 | grep TT
164 
165 //==========================================================
166 //==========================================================
167 void
168 Vertex3D::study(TVector3 r, int eveID){
169  float Z0=r.z();
170 
171  hA[0]->Fill(10);
172  hA[1]->Fill(track.size());
173  dumpPrimTracks4beamLine(Z0, eveID);
174  trackChi2QA(Z0);
175  // doExtrapolation(); /* for cov-matrix study, activate by hand */
176  if (track.size()<cut_numTrack) return;
177  hA[0]->Fill(11);
178 
179  // fill event plots .........
180  if(nHE<mxHE) {
181  TList *Lyx= hYX[nHE]->GetListOfFunctions();
182  TList *Lyz= hYZ[nHE]->GetListOfFunctions();
183  TArrow *ar=0;
184  //printf("Z0=%f, eveID=%d\n",Z0,eveID);
185  float D=2.5; // cm half length of tracklet
186  float mxPt=0;
187  for(unsigned int i=0;i<track.size();i++) {
188  DcaTrack tr=track[i]->dcaTrack;
189  float sig=tr.sigYloc;
190  float x=tr.R.x();
191  float y=tr.R.y();
192  float z=tr.R.z();
193  float px=tr.gP.x();
194  float py=tr.gP.y();
195  float pz=tr.gP.z();
196  float pt=tr.gP.Pt();
197  float pyz=sqrt(py*py+pz*pz);
198  float ux = D*px/pt; //formula uses unit vector in Y_Z plane
199  float uy = D*py/pt;
200  float vz = D*pz/pyz;//... in Y-Z plane
201  float vy = D*py/pyz;
202  float head=pt/100.;
203  int width=int(3.* (0.2*0.2)/sig/sig);
204  if(mxPt<pt) mxPt=pt;
205  // printf("tr PT=%f x=%f dx=%f y=%f dy=%f z=%f vz=%f head=%f width=%d sig=%f\n",pt, x,ux,y,uy,z,vz,head, width,sig);
206 
207  ar=new TArrow(x-ux,y-uy,x+ux,y+uy,head,"|>");
208  ar->SetLineColor(kBlue); ar->SetLineWidth(width); Lyx->Add(ar);
209 
210  ar=new TArrow(z-vz,y-vy,z+vz,y+vy,head,"|>");
211  ar->SetLineColor(kMagenta);ar->SetLineWidth(width); Lyz->Add(ar);
212  }
213 
214  // draw reco vertex 3D , dashed
215  TEllipse *el=0;
216  float Rel=0.3;
217  el=new TEllipse(r.x(),r.y(),Rel,Rel); el->SetLineColor(kRed);el->SetLineStyle(2);
218  Lyx->Add(el);
219  el=new TEllipse(r.z(),r.y(),Rel,Rel); el->SetLineColor(kBlack);el->SetLineStyle(2);
220  Lyz->Add(el);
221 
222  // try to add geant vertex to the plot
223  St_DataSet *gds=StMaker::GetChain()->GetDataSet("geant");
224  if(gds) {
225  St_g2t_vertex *g2t_ver=( St_g2t_vertex *)gds->Find("g2t_vertex");
226  if( g2t_ver) {
227  g2t_vertex_st *GVER= g2t_ver->GetTable();
228  float *GV=GVER->ge_x;
229  printf("#GGVER x=%.1f y=%.1f z=%.1f ",GV[0],GV[1],GV[2]);
230  TEllipse *el=0;
231  float Rel=0.25;
232  el=new TEllipse(GV[0],GV[1],Rel,Rel); el->SetLineColor(kRed);
233  Lyx->Add(el);
234 
235  el=new TEllipse(GV[2],GV[1],Rel,Rel); el->SetLineColor(kBlack);
236  Lyz->Add(el);
237  }
238  }
239 
240  // change title of Y_X
241  char tt[1000]; sprintf(tt,"%s, mxPt=%.1f GeV/c",hYX[nHE]->GetTitle(), mxPt);
242  hYX[nHE]->SetTitle(tt);
243 
244  // change title of Y_Z
245  sprintf(tt,"%s, eveID=%d",hYZ[nHE]->GetTitle(), eveID);
246  hYZ[nHE]->SetTitle(tt);
247 
248  // change Z_range for Y_Z
249  hYZ[nHE]->GetXaxis()->Set(2,Z0-3,Z0+3);
250  nHE++;
251  }
252 }
253 
254 
255 //==========================================================
256 void Vertex3D::dumpPrimTracks4beamLine(float z0, int eveID)
257 {
258 if (z0 || eveID) {}
259 #if 0
260  for(unsigned int i=0;i<track.size();i++) {
261  DcaTrack tr=track[i]->dcaTrack;
262  StiNodeErrs *er=&(tr.fitErr);
263  float x=tr.R.x();
264  float y=tr.R.y();
265  float z=tr.R.z();
266  float px=tr.gP.x();
267  float py=tr.gP.y();
268  float pz=tr.gP.z();
269  printf("track4beamLine %f %f %f %f %f %f %f %f %f %d %f %.1f %d \n",x,y,z,px,py,pz,er->_cYY,er->_cZY,er->_cZZ , tr.nFitPoint,tr.gChi2,z0,eveID);
270  }
271 #endif
272 }
273 
274 
275 //==========================================================
276 void
277 Vertex3D::trackChi2QA(float z0) {
278  for(unsigned int i=0;i<track.size();i++) {
279  DcaTrack tr=track[i]->dcaTrack;
280  float chi2dof=tr.gChi2;
281  hA[10]->Fill(chi2dof);
282  hA[11]->Fill(chi2dof,1./tr.gP.Pt());
283  hA[12]->Fill(chi2dof,tr.gP.Eta());
284  hA[13]->Fill(chi2dof,z0);
285  hA[14]->Fill(chi2dof,tr.nFitPoint);
286  hA[15]->Fill(chi2dof,tr.gP.Phi());
287  }
288 }
289 
290 
291 //==========================================================
292 //==========================================================
293 void
294 Vertex3D:: initHisto (TObjArray*HList){
295  TList *L=0; TLine *ln=0; float yMax=1e6; TH1* h=0;
296 
297  memset(hA,0,sizeof(hA));
298 
299  assert(HList);
300  hA[0]=new TH1F("v3D_myStat","3DVF my statistics;cases",15,0.5,15.5); // count cases
301  hA[1]=h=new TH1F("v3D_inTr","# of input tracks, 3DVF; # of tracks",30,-0.5,29.5);
302  ln=new TLine(cut_numTrack-0.5,0,cut_numTrack-0.5,yMax);
303  ln->SetLineColor(kRed); ln->SetLineStyle(2);
304  L= h->GetListOfFunctions(); L->Add(ln);
305 
306  hA[2]=h=new TH1F("v3D_inPt","PT , input tracks 3DVF; global PT (GeV) ",100,0,25);
307  L= h->GetListOfFunctions();
308  ln=new TLine(cut_pT1,0,cut_pT1,yMax);
309  ln->SetLineColor(kRed); ln->SetLineStyle(2); L->Add(ln);
310  ln=new TLine(cut_pT2,0,cut_pT2,yMax);
311  ln->SetLineColor(kRed); ln->SetLineStyle(2); L->Add(ln);
312 
313  hA[3]=h=new TH1F("v3D_inSigY","sig(Yloc) @DCA , input tracks 3DVF; sig(Ylocal) (cm) ",30,0,1.2);
314  ln=new TLine(cut_sigY,0,cut_sigY,yMax);
315  ln->SetLineColor(kRed); ln->SetLineStyle(2);
316  L= h->GetListOfFunctions(); L->Add(ln);
317 
318 
319  hA[4]=new TH1F("v3D_inSigZ","sig(Z) @DCA , input tracks 3DVF; sig(Z) (cm) ",30,0,1.5);
320 
321  hA[5]=new TH1F("v3D_myStat2","3DVF track extrapolaton statistics;cases",15,0.5,15.5); // count cases
322  // free:6-9
323 
324  hA[10]=new TH1F("v3D_chi1","#chi2/DOF best prim track candidates; #chi^{2}/DOF", 100,0,5);
325  hA[11]=new TH2F("v3D_chi2","#chi2(1/PT) best prim track candidates; #chi^{2}/DOF; 1GeV/PT ", 30,0,3,30,0,2);
326  hA[12]=new TH2F("v3D_chi3","#chi2(eta) best prim track candidates; #chi^{2}/DOF; eta ", 30,0,3,30,-2,2);
327  hA[13]=new TH2F("v3D_chi4","#chi2(zVert) best prim track candidates; #chi^{2}/DOF; zVertex(cm) ", 30,0,3,30,-150,150);
328  hA[14]=new TH2F("v3D_chi5","#chi2(nFitPoints) best prim track candidates; #chi^{2}/DOF; nFitPoints ", 30,0,3,25,0,50);
329  hA[15]=new TH2F("v3D_chi6","#chi2(phi) best prim track candidates; #chi^{2}/DOF; #phi (rad) ", 30,0,3,30,-3.2,3.2);
330 
331  for(int i=0;i<mxHA;i++)
332  if(hA[i]) HList->Add(hA[i]);
333 
334 
335  // initialize event histograms..............
336  memset(hYX,0,sizeof(hYX));
337  memset(hYZ,0,sizeof(hYZ));
338 
339  nHE=0; // clear counter
340  float mxR1=3; // cm
341  for(int i=0;i<mxHE;i++) {
342  char tt1[100], tt2[1000];
343  sprintf(tt1,"v3D_eve%dYX",i);
344  sprintf(tt2,"Y-X glob tracks at DCA, ieve=%d; X (cm); Y(cm); ",i);
345  hYX[i]=new TH2F(tt1,tt2,2,-mxR1,mxR1,2,-mxR1,mxR1);
346  HList->Add(hYX[i]);
347 
348  sprintf(tt1,"v3D_eve%dYZ",i);
349  sprintf(tt2,"Y-Z glob tracks at DCA, ieve=%d; Z (cm); Y(cm); ",i);
350  hYZ[i]=new TH2F(tt1,tt2,2,-mxR1,mxR1,2,-mxR1,mxR1);
351  HList->Add(hYZ[i]);
352 
353  }
354 }
355 }// end namespace StEvPPV
356 
StiKalmanTrackNode::propagate
int propagate(StiKalmanTrackNode *p, const StiDetector *tDet, int dir)
Propagates a track encapsulated by the given node &quot;p&quot; to the given detector &quot;tDet&quot;.
Definition: StiKalmanTrackNode.cxx:1079
StiKalmanTrack
Definition of Kalman Track.
Definition: StiKalmanTrack.h:82
DcaTrack
approximtion of track as stright line @ DCA to beamLine=0,0
Definition: TrackData.h:18
StiHit
Definition: StiHit.h:51
StiKalmanTrack::getOuterMostNode
StiKalmanTrackNode * getOuterMostNode(int qua=0) const
Accessor method returns the outer most node associated with the track.
Definition: StiKalmanTrack.h:580
StiKalmanTrack::getFitPointCount
int getFitPointCount(int detectorId=0) const
Returns the number of hits associated and used in the fit of this track.
Definition: StiKalmanTrack.cxx:970
StiKalmanTrack::getInnerMostNode
StiKalmanTrackNode * getInnerMostNode(int qua=0) const
Accessor method returns the inner most node associated with the track.
Definition: StiKalmanTrack.h:594
StiKalmanTrack::getCurvature
double getCurvature() const
Return the curvature of the track at its inner most point.
Definition: StiKalmanTrack.h:411
StiKalmanTrack::getPseudoRapidity
double getPseudoRapidity() const
Return the pseudorapidity of the track.
Definition: StiKalmanTrack.h:450
StiKalmanTrack::getChi2
double getChi2() const
Definition: StiKalmanTrack.cxx:820
StiKalmanTrack::getPhi
double getPhi() const
Return azimuthal angle at inner most point of the track.
Definition: StiKalmanTrack.h:464
TDataSet::Find
virtual TDataSet * Find(const char *path) const
Definition: TDataSet.cxx:362
StiKalmanTrack::getPt
double getPt() const
Calculates and returns the transverse momentum of the track at the inner most node.
Definition: StiKalmanTrack.h:401
Generated by   doxygen 1.8.5