StRoot: StGenericVertexMaker/StiPPVertex/Vertex3D.cxx Source File

00001 #include <string.h>
00002 #include <assert.h>
00003 #include <stdio.h>
00004 #include <StMessMgr.h>
00005 
00006 #include <TH1.h>
00007 #include <TH2.h>
00008 #include <TObjArray.h>
00009 #include <TList.h>
00010 #include <TLine.h>
00011 #include <TArrow.h>
00012 #include <TEllipse.h>
00013 
00014 #include <Sti/StiKalmanTrack.h>
00015 #include <StiUtilities/StiDebug.h>
00016 //#include <Sti/StiKalmanTrackNode.h>
00017 //#include <Sti/StiTrackContainer.h>
00018 
00019 #include <StBFChain/StBFChain.h> // for geant vertex
00020 #include <tables/St_g2t_vertex_Table.h> //for geant vertex
00021 
00022 #include "Vertex3D.h"
00023 
00024 
00025 //==========================================================
00026 //==========================================================
00027 Vertex3D::Vertex3D() {
00028 }
00029 
00030 //==========================================================
00031 //==========================================================
00032 void
00033 Vertex3D::initRun(){
00034   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;
00035  }
00036 
00037 //==========================================================
00038 //==========================================================
00039 Vertex3D::~Vertex3D(){
00040 }
00041 
00042 //==========================================================
00043 //==========================================================
00044 void
00045 Vertex3D::clearEvent(){
00046   hA[0]->Fill(1);
00047   isFound=false;
00048   track.clear();
00049   //  printf("ver3DclearEvt\n");
00050 }
00051 
00052 
00053 //==========================================================
00054 //==========================================================
00055 void
00056 Vertex3D::clearTracks(){
00057   hA[0]->Fill(2);
00058   track.clear();
00059   // printf("ver3DclearTrk\n");
00060 }
00061 
00062 
00063 //==========================================================
00064 //==========================================================
00065 void
00066 Vertex3D::addTrack(TrackData* trk){
00067   if(isValid()) return; // block adding new track after 1st vertex is found
00068   if(  !trk->mTpc && !trk->mBtof && !trk->mCtb && !trk->mBemc && trk->mEemc ) return;  // dongx
00069   hA[0]->Fill(3);
00070   // track is matched to any fast detector
00071 
00072   float pt=trk->dcaTrack.gP.Pt();
00073   hA[2]->Fill(pt);
00074   if (pt <cut_pT1) return;
00075   hA[0]->Fill(4);
00076   if (pt >cut_pT2) return;
00077   hA[0]->Fill(5);
00078 
00079   float sY=trk->dcaTrack.sigYloc; //track extrapolation error in transverse direction
00080   hA[3]->Fill(sY);
00081   if(sY>cut_sigY) return;
00082   hA[0]->Fill(6);
00083 
00084   hA[4]->Fill(trk->dcaTrack.sigZ); 
00085   track.push_back(trk); 
00086 }
00087 
00088 //==========================================================
00089 //==========================================================
00090 void
00091 Vertex3D::doExtrapolation(){ // study track cov matrix for individual tracks using my macro plTrCov.C 
00092   for(uint i=0;i<track.size();i++) {
00093      const StiKalmanTrack* tr=track[i]->mother;
00094      hA[5]->Fill(1);
00095      if( fabs(tr->getChi2()-1.)>0.8) continue;
00096      hA[5]->Fill(2);
00097      if( tr->getFitPointCount()<20) continue;
00098      hA[5]->Fill(3);
00099      if( fabs(tr->getPseudoRapidity())>0.5) continue;
00100      hA[5]->Fill(4);
00101      if( fabs(tr->getPt())<1.) continue;
00102      hA[5]->Fill(5);
00103       
00104      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());
00105      printf("T: track position along nodes, N=%d\n",tr->getFitPointCount());
00106      StiKTNIterator tNode = tr->rbegin();
00107      StiKTNIterator eNode = tr->rend();
00108      int iN=0;
00109      float minLx=9999, maxLx=0.; // local coordinates of track
00110      for (;tNode!=eNode;++tNode){
00111        StiKalmanTrackNode *node = &(*tNode);
00112        if(!node->isValid())      continue;
00113      
00114        StiHit *stiHit = node->getHit();
00115        if (!stiHit)              continue;
00116          
00117        if (node->getChi2()>1000) continue;
00118        if (!node->isFitted())    continue;
00119        iN++;
00120        float Lx=node->x();
00121        float Ly=node->y();
00122        float Lz=node->z();
00123        float sLy=sqrt(node->fitErrs()._cYY);
00124        float sLz=sqrt(node->fitErrs()._cZZ);
00125        float Rxy=sqrt(Lx*Lx+Ly*Ly);
00126        if(minLx>Lx) minLx=Lx;
00127        if(maxLx<Lx) maxLx=Lx;
00128        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());
00129          printf("TTm %.2f    %.2f %.3f   %.2f %.3f  %f\n",node->x_g(),node->y_g(),sLy,node->z_g(),sLz,node->getAlpha());
00130 
00131      } // end of nodes w/ hits
00132  
00133      for(int iDir=0;iDir<2;iDir++) {
00134        if(iDir==0)  printf("extrapolate track toward vertex, innLx=%.1f\n",minLx);
00135        if(iDir==1)  printf("extrapolate track toward BSMD, outLx=%.1f\n",maxLx);       
00136        for(int iI=0;iI<17;iI++) {
00137          float Lx=100.;StiKalmanTrackNode* node=0; int ret=-2;
00138          StiKalmanTrack track2=*tr; // clone track to allow extrapolation
00139          if(iDir==0) {// toward vertex
00140            Lx=minLx-iI*5;
00141            node=track2.getInnerMostNode();
00142            ret=node->propagate(Lx,kCylindrical,kOutsideIn);
00143          } else {
00144            Lx=maxLx+iI*2; if(Lx>250) break;
00145            node=track2.getOuterMostNode();
00146            ret=node->propagate(Lx,kCylindrical,kInsideOut);
00147          }
00148          if(iI==0)cout<<"TT track extrap node:"<< *node;
00149          if(ret<0) {printf("prop in x=%f ret=%d, skip\n",Lx,ret); continue;}
00150          float Ly=node->y();
00151          float Lz=node->z();
00152          if(fabs(Ly)>40.) break; // to not cross sector boundary ??
00153          if(iDir==1 && fabs(Ly)>36.) break; // to not cross sector boundary ??
00154          node->propagateError();
00155          float sLy=sqrt(node->fitErrs()._cYY);
00156          float sLz=sqrt(node->fitErrs()._cZZ);
00157          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());
00158          printf("TT%d %.2f   %.2f %.3f   %.2f %.3f  %f\n",iDir,node->x_g(),node->y_g(),sLy,node->z_g(),sLz,node->getAlpha());
00159        }
00160      }
00161      
00162   }  
00163 }
00164 //head -n 4224 Lda83-618 | tail -n 150 | grep TT
00165 
00166 //==========================================================
00167 //==========================================================
00168 void
00169 Vertex3D::study(TVector3 r, int eveID){
00170   float Z0=r.z();
00171   
00172   hA[0]->Fill(10);
00173   hA[1]->Fill(track.size());
00174   dumpPrimTracks4beamLine(Z0, eveID);
00175   trackChi2QA(Z0);
00176   // doExtrapolation(); /* for cov-matrix study, activate by hand */
00177   if (track.size()<cut_numTrack)  return;
00178   hA[0]->Fill(11);
00179 
00180   // fill event plots .........
00181   if(nHE<mxHE) {
00182     TList *Lyx= hYX[nHE]->GetListOfFunctions(); 
00183     TList *Lyz= hYZ[nHE]->GetListOfFunctions(); 
00184     TArrow *ar=0;
00185     //printf("Z0=%f, eveID=%d\n",Z0,eveID);
00186     float D=2.5; // cm half length of tracklet
00187     float mxPt=0;
00188     for(uint i=0;i<track.size();i++) {
00189       DcaTrack  tr=track[i]->dcaTrack;
00190       float sig=tr.sigYloc;
00191       float x=tr.R.x();
00192       float y=tr.R.y();
00193       float z=tr.R.z();
00194       float px=tr.gP.x();
00195       float py=tr.gP.y();
00196       float pz=tr.gP.z();
00197       float pt=tr.gP.Pt();
00198       float pyz=sqrt(py*py+pz*pz);
00199       float ux = D*px/pt; //formula uses unit vector in Y_Z plane
00200       float uy = D*py/pt;
00201       float vz = D*pz/pyz;//... in Y-Z plane
00202       float vy = D*py/pyz;
00203       float head=pt/100.;
00204       int width=int(3.* (0.2*0.2)/sig/sig);
00205       if(mxPt<pt) mxPt=pt;
00206       // 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);
00207 
00208       ar=new TArrow(x-ux,y-uy,x+ux,y+uy,head,"|>");  
00209       ar->SetLineColor(kBlue); ar->SetLineWidth(width); Lyx->Add(ar);
00210 
00211       ar=new TArrow(z-vz,y-vy,z+vz,y+vy,head,"|>"); 
00212       ar->SetLineColor(kMagenta);ar->SetLineWidth(width); Lyz->Add(ar);
00213     }
00214 
00215     // draw reco vertex 3D , dashed
00216     TEllipse *el=0;
00217     float Rel=0.3;
00218     el=new TEllipse(r.x(),r.y(),Rel,Rel); el->SetLineColor(kRed);el->SetLineStyle(2);
00219     Lyx->Add(el);
00220     el=new TEllipse(r.z(),r.y(),Rel,Rel); el->SetLineColor(kBlack);el->SetLineStyle(2);
00221     Lyz->Add(el);
00222     
00223     // try to add geant vertex to the plot
00224     St_DataSet *gds=StMaker::GetChain()->GetDataSet("geant");
00225     if(gds) {
00226       St_g2t_vertex  *g2t_ver=( St_g2t_vertex *)gds->Find("g2t_vertex");
00227       if( g2t_ver) {
00228         g2t_vertex_st *GVER= g2t_ver->GetTable();
00229         float *GV=GVER->ge_x;
00230         printf("#GGVER  x=%.1f y=%.1f z=%.1f ",GV[0],GV[1],GV[2]); 
00231         TEllipse *el=0;
00232         float Rel=0.25;
00233         el=new TEllipse(GV[0],GV[1],Rel,Rel); el->SetLineColor(kRed);
00234         Lyx->Add(el);
00235         
00236         el=new TEllipse(GV[2],GV[1],Rel,Rel); el->SetLineColor(kBlack);
00237         Lyz->Add(el);
00238       }
00239     }
00240     
00241     // change title of Y_X
00242     char tt[1000]; sprintf(tt,"%s, mxPt=%.1f GeV/c",hYX[nHE]->GetTitle(), mxPt);
00243     hYX[nHE]->SetTitle(tt);
00244     
00245     // change title of Y_Z
00246     sprintf(tt,"%s, eveID=%d",hYZ[nHE]->GetTitle(), eveID);
00247     hYZ[nHE]->SetTitle(tt);
00248 
00249     // change Z_range for Y_Z
00250     hYZ[nHE]->GetXaxis()->Set(2,Z0-3,Z0+3);
00251     nHE++;
00252   }
00253 }
00254 
00255 
00256 //==========================================================
00257 void
00258 Vertex3D::dumpPrimTracks4beamLine(float z0, int eveID) {
00259   for(uint i=0;i<track.size();i++) {
00260       DcaTrack  tr=track[i]->dcaTrack;
00261       StiNodeErrs *er=&(tr.fitErr);
00262       float x=tr.R.x();
00263       float y=tr.R.y();
00264       float z=tr.R.z();
00265       float px=tr.gP.x();
00266       float py=tr.gP.y();
00267       float pz=tr.gP.z();
00268       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);
00269   }
00270 }
00271 
00272 
00273 //==========================================================
00274 void
00275 Vertex3D::trackChi2QA(float z0) {
00276   for(uint i=0;i<track.size();i++) {
00277       DcaTrack  tr=track[i]->dcaTrack;
00278        float chi2dof=tr.gChi2;
00279       hA[10]->Fill(chi2dof);
00280       hA[11]->Fill(chi2dof,1./tr.gP.Pt());
00281       hA[12]->Fill(chi2dof,tr.gP.Eta());
00282       hA[13]->Fill(chi2dof,z0);
00283       hA[14]->Fill(chi2dof,tr.nFitPoint);
00284       hA[15]->Fill(chi2dof,tr.gP.Phi());
00285   }
00286 }
00287 
00288 
00289 //==========================================================
00290 //==========================================================
00291 void
00292 Vertex3D:: initHisto (TObjArray*HList){
00293   TList *L=0;  TLine *ln=0; float yMax=1e6; TH1* h=0;
00294   
00295   memset(hA,0,sizeof(hA));
00296  
00297   assert(HList);
00298   hA[0]=new TH1F("v3D_myStat","3DVF my statistics;cases",15,0.5,15.5); // count cases
00299   hA[1]=h=new TH1F("v3D_inTr","# of input tracks, 3DVF; # of tracks",30,-0.5,29.5);
00300   ln=new TLine(cut_numTrack-0.5,0,cut_numTrack-0.5,yMax);
00301   ln->SetLineColor(kRed); ln->SetLineStyle(2);
00302   L= h->GetListOfFunctions();  L->Add(ln);
00303 
00304   hA[2]=h=new TH1F("v3D_inPt","PT , input tracks  3DVF; global PT (GeV) ",100,0,25);
00305   L= h->GetListOfFunctions(); 
00306   ln=new TLine(cut_pT1,0,cut_pT1,yMax);
00307   ln->SetLineColor(kRed); ln->SetLineStyle(2); L->Add(ln);
00308   ln=new TLine(cut_pT2,0,cut_pT2,yMax);
00309   ln->SetLineColor(kRed); ln->SetLineStyle(2); L->Add(ln);
00310 
00311   hA[3]=h=new TH1F("v3D_inSigY","sig(Yloc) @DCA , input tracks  3DVF; sig(Ylocal) (cm) ",30,0,1.2);
00312   ln=new TLine(cut_sigY,0,cut_sigY,yMax);
00313   ln->SetLineColor(kRed); ln->SetLineStyle(2);
00314   L= h->GetListOfFunctions();  L->Add(ln);
00315 
00316 
00317   hA[4]=new TH1F("v3D_inSigZ","sig(Z) @DCA , input tracks  3DVF; sig(Z) (cm) ",30,0,1.5);
00318 
00319   hA[5]=new TH1F("v3D_myStat2","3DVF track extrapolaton statistics;cases",15,0.5,15.5); // count cases
00320   // free:6-9
00321 
00322   hA[10]=new TH1F("v3D_chi1","#chi2/DOF best prim track candidates; #chi^{2}/DOF", 100,0,5);
00323   hA[11]=new TH2F("v3D_chi2","#chi2(1/PT) best prim track candidates; #chi^{2}/DOF; 1GeV/PT ", 30,0,3,30,0,2);
00324   hA[12]=new TH2F("v3D_chi3","#chi2(eta) best prim track candidates; #chi^{2}/DOF; eta ", 30,0,3,30,-2,2);
00325   hA[13]=new TH2F("v3D_chi4","#chi2(zVert) best prim track candidates; #chi^{2}/DOF; zVertex(cm) ", 30,0,3,30,-150,150);  
00326   hA[14]=new TH2F("v3D_chi5","#chi2(nFitPoints) best prim track candidates; #chi^{2}/DOF; nFitPoints ", 30,0,3,25,0,50);
00327   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);
00328 
00329   for(int i=0;i<mxHA;i++) 
00330     if(hA[i]) HList->Add(hA[i]);
00331 
00332 
00333   // initialize event histograms..............
00334   memset(hYX,0,sizeof(hYX));
00335   memset(hYZ,0,sizeof(hYZ));
00336 
00337   nHE=0; // clear counter
00338   float mxR1=3; // cm
00339   for(int i=0;i<mxHE;i++) {
00340     char tt1[100], tt2[1000];
00341     sprintf(tt1,"v3D_eve%dYX",i);
00342     sprintf(tt2,"Y-X glob tracks at DCA, ieve=%d; X (cm); Y(cm); ",i);
00343     hYX[i]=new TH2F(tt1,tt2,2,-mxR1,mxR1,2,-mxR1,mxR1);
00344     HList->Add(hYX[i]);
00345 
00346     sprintf(tt1,"v3D_eve%dYZ",i);
00347     sprintf(tt2,"Y-Z glob tracks at DCA, ieve=%d; Z (cm); Y(cm); ",i);
00348     hYZ[i]=new TH2F(tt1,tt2,2,-mxR1,mxR1,2,-mxR1,mxR1);
00349     HList->Add(hYZ[i]);
00350 
00351   }
00352 }
00353