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 <St_base/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 <Sti/StiKalmanTrack.h>
14 #include <StiUtilities/StiDebug.h>
15 //#include <Sti/StiKalmanTrackNode.h>
16 //#include <Sti/StiTrackContainer.h>
17 
18 #include <StBFChain/StBFChain.h> // for geant vertex
19 #include <tables/St_g2t_vertex_Table.h> //for geant vertex
20 
21 #include "StGenericVertexMaker/StiPPVertex/Vertex3D.h"
22 
23 
24 //==========================================================
25 //==========================================================
26 Vertex3D::Vertex3D() {
27 }
28 
29 //==========================================================
30 //==========================================================
31 void
32 Vertex3D::initRun(){
33  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;
34  }
35 
36 //==========================================================
37 //==========================================================
38 Vertex3D::~Vertex3D(){
39 }
40 
41 //==========================================================
42 //==========================================================
43 void
44 Vertex3D::clearEvent(){
45  hA[0]->Fill(1);
46  isFound=false;
47  track.clear();
48  // printf("ver3DclearEvt\n");
49 }
50 
51 
52 //==========================================================
53 //==========================================================
54 void
55 Vertex3D::clearTracks(){
56  hA[0]->Fill(2);
57  track.clear();
58  // printf("ver3DclearTrk\n");
59 }
60 
61 
62 //==========================================================
63 //==========================================================
64 void
65 Vertex3D::addTrack(TrackData* trk){
66  if(isValid()) return; // block adding new track after 1st vertex is found
67  if( !trk->mTpc && !trk->mBtof && !trk->mCtb && !trk->mBemc && trk->mEemc ) return;
68  hA[0]->Fill(3);
69  // track is matched to any fast detector
70 
71  float pt=trk->dcaTrack.gP.Pt();
72  hA[2]->Fill(pt);
73  if (pt <cut_pT1) return;
74  hA[0]->Fill(4);
75  if (pt >cut_pT2) return;
76  hA[0]->Fill(5);
77 
78  float sY=trk->dcaTrack.sigYloc; //track extrapolation error in transverse direction
79  hA[3]->Fill(sY);
80  if(sY>cut_sigY) return;
81  hA[0]->Fill(6);
82 
83  hA[4]->Fill(trk->dcaTrack.sigZ);
84  track.push_back(trk);
85 }
86 
87 //==========================================================
88 //==========================================================
89 void
90 Vertex3D::doExtrapolation(){ // study track cov matrix for individual tracks using my macro plTrCov.C
91  for(unsigned int i=0;i<track.size();i++) {
92  const StiKalmanTrack* tr=track[i]->getMother<StiKalmanTrack>();
93  hA[5]->Fill(1);
94  if( fabs(tr->getChi2()-1.)>0.8) continue;
95  hA[5]->Fill(2);
96  if( tr->getFitPointCount()<20) continue;
97  hA[5]->Fill(3);
98  if( fabs(tr->getPseudoRapidity())>0.5) continue;
99  hA[5]->Fill(4);
100  if( fabs(tr->getPt())<1.) continue;
101  hA[5]->Fill(5);
102 
103  const StiKalmanTrackNode* node = tr->getInnOutMostNode(0, 3);
104  StThreeVector<double> trackNodeCoord(node->x_g(),node->y_g(),node->z_g());
105 
106  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(), trackNodeCoord.perp(),tr->getFitPointCount(),tr->getPseudoRapidity(),tr->getChi2());
107  printf("T: track position along nodes, N=%d\n",tr->getFitPointCount());
108  StiKTNIterator tNode = tr->rbegin();
109  StiKTNIterator eNode = tr->rend();
110  int iN=0;
111  float minLx=9999, maxLx=0.; // local coordinates of track
112  for (;tNode!=eNode;++tNode){
113  StiKalmanTrackNode *node = &(*tNode);
114  if(!node->isValid()) continue;
115 
116  StiHit *stiHit = node->getHit();
117  if (!stiHit) continue;
118 
119  if (node->getChi2()>1000) continue;
120  if (!node->isFitted()) continue;
121  iN++;
122  float Lx=node->x();
123  float Ly=node->y();
124  float Lz=node->z();
125  float sLy=sqrt(node->fitErrs()._cYY);
126  float sLz=sqrt(node->fitErrs()._cZZ);
127  float Rxy=sqrt(Lx*Lx+Ly*Ly);
128  if(minLx>Lx) minLx=Lx;
129  if(maxLx<Lx) maxLx=Lx;
130  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());
131  printf("TTm %.2f %.2f %.3f %.2f %.3f %f\n",node->x_g(),node->y_g(),sLy,node->z_g(),sLz,node->getAlpha());
132 
133  } // end of nodes w/ hits
134 
135  for(int iDir=0;iDir<2;iDir++) {
136  if(iDir==0) printf("extrapolate track toward vertex, innLx=%.1f\n",minLx);
137  if(iDir==1) printf("extrapolate track toward BSMD, outLx=%.1f\n",maxLx);
138  for(int iI=0;iI<17;iI++) {
139  float Lx=100.;StiKalmanTrackNode* node=0; int ret=-2;
140  StiKalmanTrack track2=*tr; // clone track to allow extrapolation
141  if(iDir==0) {// toward vertex
142  Lx=minLx-iI*5;
143  node=track2.getInnerMostNode();
144  ret=node->propagate(Lx,kCylindrical,kOutsideIn);
145  } else {
146  Lx=maxLx+iI*2; if(Lx>250) break;
147  node=track2.getOuterMostNode();
148  ret=node->propagate(Lx,kCylindrical,kInsideOut);
149  }
150  if(iI==0)cout<<"TT track extrap node:"<< *node;
151  if(ret<0) {printf("prop in x=%f ret=%d, skip\n",Lx,ret); continue;}
152  float Ly=node->y();
153  float Lz=node->z();
154  if(fabs(Ly)>40.) break; // to not cross sector boundary ??
155  if(iDir==1 && fabs(Ly)>36.) break; // to not cross sector boundary ??
156  node->propagateError();
157  float sLy=sqrt(node->fitErrs()._cYY);
158  float sLz=sqrt(node->fitErrs()._cZZ);
159  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());
160  printf("TT%d %.2f %.2f %.3f %.2f %.3f %f\n",iDir,node->x_g(),node->y_g(),sLy,node->z_g(),sLz,node->getAlpha());
161  }
162  }
163 
164  }
165 }
166 //head -n 4224 Lda83-618 | tail -n 150 | grep TT
167 
168 //==========================================================
169 //==========================================================
170 void
171 Vertex3D::study(TVector3 r, int eveID){
172  float Z0=r.z();
173 
174  hA[0]->Fill(10);
175  hA[1]->Fill(track.size());
176  dumpPrimTracks4beamLine(Z0, eveID);
177  trackChi2QA(Z0);
178  // doExtrapolation(); /* for cov-matrix study, activate by hand */
179  if (track.size()<cut_numTrack) return;
180  hA[0]->Fill(11);
181 
182  // fill event plots .........
183  if(nHE<mxHE) {
184  TList *Lyx= hYX[nHE]->GetListOfFunctions();
185  TList *Lyz= hYZ[nHE]->GetListOfFunctions();
186  TArrow *ar=0;
187  //printf("Z0=%f, eveID=%d\n",Z0,eveID);
188  float D=2.5; // cm half length of tracklet
189  float mxPt=0;
190  for(unsigned int i=0;i<track.size();i++) {
191  DcaTrack tr=track[i]->dcaTrack;
192  float sig=tr.sigYloc;
193  float x=tr.R.x();
194  float y=tr.R.y();
195  float z=tr.R.z();
196  float px=tr.gP.x();
197  float py=tr.gP.y();
198  float pz=tr.gP.z();
199  float pt=tr.gP.Pt();
200  float pyz=sqrt(py*py+pz*pz);
201  float ux = D*px/pt; //formula uses unit vector in Y_Z plane
202  float uy = D*py/pt;
203  float vz = D*pz/pyz;//... in Y-Z plane
204  float vy = D*py/pyz;
205  float head=pt/100.;
206  int width=int(3.* (0.2*0.2)/sig/sig);
207  if(mxPt<pt) mxPt=pt;
208  // 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);
209 
210  ar=new TArrow(x-ux,y-uy,x+ux,y+uy,head,"|>");
211  ar->SetLineColor(kBlue); ar->SetLineWidth(width); Lyx->Add(ar);
212 
213  ar=new TArrow(z-vz,y-vy,z+vz,y+vy,head,"|>");
214  ar->SetLineColor(kMagenta);ar->SetLineWidth(width); Lyz->Add(ar);
215  }
216 
217  // draw reco vertex 3D , dashed
218  TEllipse *el=0;
219  float Rel=0.3;
220  el=new TEllipse(r.x(),r.y(),Rel,Rel); el->SetLineColor(kRed);el->SetLineStyle(2);
221  Lyx->Add(el);
222  el=new TEllipse(r.z(),r.y(),Rel,Rel); el->SetLineColor(kBlack);el->SetLineStyle(2);
223  Lyz->Add(el);
224 
225  // try to add geant vertex to the plot
226  St_DataSet *gds=StMaker::GetChain()->GetDataSet("geant");
227  if(gds) {
228  St_g2t_vertex *g2t_ver=( St_g2t_vertex *)gds->Find("g2t_vertex");
229  if( g2t_ver) {
230  g2t_vertex_st *GVER= g2t_ver->GetTable();
231  float *GV=GVER->ge_x;
232  printf("#GGVER x=%.1f y=%.1f z=%.1f ",GV[0],GV[1],GV[2]);
233  TEllipse *el=0;
234  float Rel=0.25;
235  el=new TEllipse(GV[0],GV[1],Rel,Rel); el->SetLineColor(kRed);
236  Lyx->Add(el);
237 
238  el=new TEllipse(GV[2],GV[1],Rel,Rel); el->SetLineColor(kBlack);
239  Lyz->Add(el);
240  }
241  }
242 
243  // change title of Y_X
244  char tt[1000]; sprintf(tt,"%s, mxPt=%.1f GeV/c",hYX[nHE]->GetTitle(), mxPt);
245  hYX[nHE]->SetTitle(tt);
246 
247  // change title of Y_Z
248  sprintf(tt,"%s, eveID=%d",hYZ[nHE]->GetTitle(), eveID);
249  hYZ[nHE]->SetTitle(tt);
250 
251  // change Z_range for Y_Z
252  hYZ[nHE]->GetXaxis()->Set(2,Z0-3,Z0+3);
253  nHE++;
254  }
255 }
256 
257 
258 //==========================================================
259 void
260 Vertex3D::dumpPrimTracks4beamLine(float z0, int eveID) {
261  for(unsigned int i=0;i<track.size();i++) {
262  DcaTrack tr=track[i]->dcaTrack;
263  StiNodeErrs *er=&(tr.fitErr);
264  float x=tr.R.x();
265  float y=tr.R.y();
266  float z=tr.R.z();
267  float px=tr.gP.x();
268  float py=tr.gP.y();
269  float pz=tr.gP.z();
270  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);
271  }
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 
StiKalmanTrack::getInnOutMostNode
StiKalmanTrackNode * getInnOutMostNode(int inot, int qua) const
Same for getNNodes(qua)
Definition: StiKalmanTrack.cxx:1195
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.h
Definition of Kalman Track.
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