/*:>-------------------------------------------------------------------- **: FILE: mMvdPseudoTrk.cc **: HISTORY: **: 00jan93-v000a-hpl- Created by stic Version **: Id: idl.y,v 1.17 1997/03/25 19:22:52 ward Exp **:<------------------------------------------------------------------*/ #include "mMvdPseudoTrk.h" #include #include #include #include "PhHistogramFactory.hh" #include void find_closest_center(float *data,int Ll,float *max,int *channel); void refine_resolution(float *data,int Ll,float step, float z1, int ichmax,float *sumx,float *sumy); void proj_histogram(DMVDGEO_ST *Geo, TABLE_HEAD_ST *Clmp_h,DMVDCLMP_ST *Clmp , float *data, float x1, float step, int ndim, float z_vtx,int row); void find_range(float z_pos,float z_vtx,float delz,DMVDGEO_ST *Geo, float *lower,float *upper); #define ZRESOLUTION 0.030 /* binsize in Z for vertex search histogram */ #define NUMDATA 100000 /* max number of bins in array used in z vertex */ /* finding the real number should be more like 6000 */ #define SOFTWARE_VERSION 1.01 #define LOUD 3 #define VERY_LOUD 4 long mMvdPseudoTrk_( TABLE_HEAD_ST *Geo_h, DMVDGEO_ST *Geo , TABLE_HEAD_ST *Par_h, DMVDPAR_ST *Par , TABLE_HEAD_ST *bPar_h, DMVBPAR_ST *bPar , TABLE_HEAD_ST *Nclmp_h, DMVDNCLMP_ST *Nclmp , TABLE_HEAD_ST *Clmp_h, DMVDCLMP_ST *Clmp , TABLE_HEAD_ST *Cntrl_h, DMVDTRIGCONTROL_ST *Cntrl , TABLE_HEAD_ST *pPar_h, DMVDPSEUDOPAR_ST *pPar , TABLE_HEAD_ST *VertexOut_h, DMVDVERTEXOUT_ST *VertexOut , TABLE_HEAD_ST *Io_h, DMVDIO_ST *Io ) { /*:>-------------------------------------------------------------------- **: ROUTINE : mMvdPseudoTrk_ **: **: DESCRIPTION: Code to find the vertex position using the **: "pseudo-tracking" method. Find the vertex location **: by projecting all possible track combinations towards **: a z position for the vertex. Then take all of these **: z positions and choose the "most popular" value. **: An array of projected track candidate combinations **: with dimension (Lz) is used. The bin-size in this **: array can potentially limit the vertex finding **: resolution. The calculated z vertex position is in **: VertexOut->vertex[2] (a staf table entry). **: **: AUTHOR : J.H.Park, S.S.Ryu (YonSei Univ.) (Nov 10 97) **: Translated from the original Fortran version VER_ZTRK.f **: written by J.P.Sullivan **: **: ARGUMENTS : **: IN: **: Geo - mvd geometry parameters from GEANT **: Geo_h - header Structure for Geo **: Par - mvd uncalibration parameters **: Par_h - header Structure for Par **: bPar - barrel detector geometry parameters **: bPar_h - header Structure for bPar **: Nclmp - # clump on each row of barrel detector **: Nclmp_h - header Structure for Nclmp **: Clmp - size and z position of each clump **: Clmp_h - header Structure for Clmp **: Io - contains verbose parameter **: Io_h - header structure for Io **: INOUT: **: VertexOut - vertex-finding output table **: VertexOut_h- header Structure for VertexOut **: pPar - pseudo tracking parameters **: pPar_h - header Structure for pPar **: OUT: **: RETURNS: STAF Condition Value **: Revision History **: Name Date Comment **: J.H.Park Feb. 10, 1998 controll execution of this module **: by the value of choice **: J.P.Sullivan Apr. 6, 1998 add code to calculate and histogram **: cpu time used. **: J.P.Sullivan Apr. 16, 1998 replace printf with Mvdmess calls, **: change default IEXPND from 3 to 2. **: J.P.Sullivan Apr. 28, 1998 remove header table from calling **: sequence and remove one debug **: message which used it and two lines **: of commented out code which used **: it (for real vertex position). **: J.P.Sullivan Jun. 15, 1998 I made this change on or about **: 26-May-98. Make the binsize easier **: to change for pseudotracking. It **: is now controlled by ZRESOLUTION. **: In the process add some comments. **: J.P.Sullivan July 21, 1998 Changed output a little, fixed **: some errors in typecasting, which **: were probably introduced in the **: June 15, 1998 changes. **: J.P.Sullivan July 29, 1998 Removed STAFCV_BAD return code **: J.P.Sullivan July 31, 1998 Removed use of hbook, vertex **: finding uses local "histogram" **: arrays x_data, y_data, z_data. **: Change to C++ (really very few **: changes were made) and replace **: mMvdMess_ calls with cout. Remove **: dMvdHst from calling sequence. **: Change histograms to use generic **: histogram interface. **: J.P.Sullivan 10-Aug-98 Add dMvdVertexOut table and code to fill it. **: Remove mMvdPseudoTrk table and the code **: which previously filled it. Use a local **: array instead of a histogram to find the **: vertex. Add code to use staf generic **: histograms for diagnostic histograms. **: Remove hbook. Replace mMvdMess with cout. **: Fixed a 1-bin bug in vertex location (i.e. **: about 300 microns). **: J.P.Sullivan 23-Sep-98 Add dMvdIo to calling sequence and use it to **: control volume of output **: S.S.Ryu 06-Jul-99 changes the variable "Cntl" to "Cntrl" to make **: this variable consistent throughout the MVD **: codes. **:>------------------------------------------------------------------*/ float xofftmp1,xofftmp2,yofftmp1,yofftmp2; float dxstep,dystep,dzstep; float zproj,ymax,sumx,sumy,zval; float ratio1,ratio2; float z2val,z1,z2,z1min,z1max; float xtmp; int jj,ii,i,j,k,ichmax,i1,i2; int ibin; unsigned innerindex=0; unsigned outerindex=0; float cputime,cpu_in,cpu_out; /* for tests of cpu time use */ static short ifirst=0; static int ncalls=0; short ErrorCode; float Xvertex, Yvertex, Zvertex; int NumOut; int Lx,Ly,Lz; int nbin; float *x_data; float *y_data; float *z_data; float x_vertex,y_vertex; float outer_radius; static PhH1 *CpuPerEvP; static PhH2 *CpuVsOccP; static PhH1 *PseudoX; static PhH1 *PseudoY; static PhH1 *PseudoZ; PhHistogramFactory *hf = PhHistogramFactory::instance(); if ( Cntrl->choice ==1 ) return STAFCV_OK; /* return if this algorithm is not selected */ cpu_in=(float)clock()/(float)CLOCKS_PER_SEC; if ( Io->verbose < LOUD ) cout << "\nmMvdPseudoTrk{" << endl; /* ***************************************************************** */ /* The first step will be to initialize (mostly zeroes) the contents */ /* of the output table. */ ncalls++; NumOut = VertexOut_h->nok; /* This saves the number of the row we */ /* are currently filling */ VertexOut_h->nok++; for ( i=0 ; i<3 ; i++ ) { VertexOut[NumOut].vertex[i]=0.; for ( j=0 ; j<3 ; j++ ) VertexOut[NumOut].vertexerr[i][j]=0.; } VertexOut[NumOut].softvertex=SOFTWARE_VERSION; VertexOut[NumOut].ErrorCode=0; /* 1 means OK, negative means error */ VertexOut[NumOut].vertexcl=1.; /* use 1 until there is a better number */ VertexOut[NumOut].id=ncalls; /* N event would be better * ErrorCode=0; /* ******************************************************************** */ if ( pPar_h->nok==0) pPar->IEXPND=2; /* Z resolution is partly controlled by the binsize */ /* expand search range by pPar->IEXPND */ z2 = pPar->IEXPND*Par->zoff; /* maximum z to look for vertex */ z1 = -z2; /* minimum z to look for vertex */ Lz = (int)((z2-z1)/ZRESOLUTION);/* N channels in histogram */ if ( Lz > NUMDATA ) Lz=NUMDATA; /* make sure it fits in the array */ z_data = (float *)calloc(Lz,sizeof(float)); /* allocate the array used*/ if ( z_data==NULL ) { if ( Io->verbose < VERY_LOUD ) cout <<"\t" << " Z vertex : out of memory ========" << endl; VertexOut[NumOut].ErrorCode=-5; return STAFCV_OK; /* should really return STAFCV_BAD */ } dzstep = (z2-z1)/(float)Lz; /* histogram binsize */ /* x,y vertex position resolution is limited from Lx,Ly,Lz values which are used as bin number */ Lx=Ly=(int)((Geo->vislr+Geo->dvisl[0])*100.); ratio1=(Geo->vislr+Geo->dvisl[0]) /((Geo->voslr+Geo->dvosl[0])-(Geo->vislr+Geo->dvisl[0])); ratio2=(Geo->vislr+Geo->dvisl[0])/(Geo->voslr+Geo->dvosl[0]); /* xofftmp1 is lower bound of x vertex position xofftmp2 is upper bound of x vertex position as yofftmp1,yofftmp2 are used in determining bounds of the y vertex position */ xofftmp1=(Geo->voslr+Geo->dvosl[0])-(Geo->vislr+Geo->dvisl[0]); xofftmp2=(Geo->voslr+Geo->dvosl[0])+(Geo->vislr+Geo->dvisl[0]); yofftmp1=xofftmp1; yofftmp2=xofftmp2; /* bin size is 200micron */ dxstep= (xofftmp2-xofftmp1)/(float)Lx; dystep= (yofftmp2-yofftmp1)/(float)Ly; /* Histogram setup */ hf->cd(); hf->cd("MVD"); if ( ifirst==0 ) { ifirst=1; CpuPerEvP = hf->CreateH1F("CpuPerEvP","CPU sec/ev (pseudotracking)", 500, 0., 5.); CpuVsOccP = hf->CreateH2F("CpuVsOccP","CPU sec/ev vs. occupancy", 20,0.,1.,100,0.,5.); PseudoX = hf->CreateH1F("PseudoX","hist for X vertex",Lx,xofftmp1,xofftmp2); PseudoY = hf->CreateH1F("PseudoY","hist for Y vertex",Ly,yofftmp1,yofftmp2); PseudoZ = hf->CreateH1F("PseudoZ","hist for Z vertex",Lz,z1,z2); } else { PseudoX->Reset(); PseudoY->Reset(); PseudoZ->Reset(); } for ( ibin=0; ibinNPHUSE;++j) { /* innerindex and outerindex point to the first clump in the current */ /* row (index for rows in j) in the inner and outer shells */ innerindex=outerindex=0; for (k=0;k< bPar->NPHUSE+j;++k) outerindex+=Nclmp[k].num; /*outer barrel Clmp index */ for (k=0;kNPHUSE+j].num >0 && Nclmp[j].num >0) { /* the following loop is over all clumps in the outer barrel */ for (jj=0;jjNPHUSE+j].num;++jj){ z2val=Clmp[outerindex+jj].zavg; /* the z of this outer shell clump */ z1min=z1+ratio2*(z2val-z1); /* min z to consider in inner shell*/ z1max=z2+ratio2*(z2val-z2); /* max z to consider in inner shell*/ if ( Clmp[outerindex+jj].shell != 1 ) { if ( Io->verbose < VERY_LOUD ) cout <<"\t" << "*** ERROR this should be shell 1 but it is shell " << Clmp[outerindex+jj].shell << endl; } /* the following loop is over all hits in the inner barrel */ /* for each hit in z the range that could possibly projected to*/ /* the selected range of vertex positions, project to the z at */ /* r=0 and put it in a histogram */ for (ii=0;iiverbose < VERY_LOUD ) cout <<"\t" << "*** ERROR shell 0 row=" << Clmp[innerindex+ii].row << " shell 1 row=" << Clmp[outerindex+jj].row << endl; if ( Clmp[innerindex+ii].shell != 0 ) if ( Io->verbose < VERY_LOUD ) cout <<"\t" << "*** ERROR this should be shell 0 but it is shell" << Clmp[innerindex+ii].shell << endl; if(zval > z1min && zval < z1max) { zproj=zval-ratio1*(z2val-zval); ibin = (int)((zproj-z1)/dzstep); if ( (0<=ibin) && (ibin0 ) { refine_resolution(z_data,Lz,dzstep,z1,ichmax,&sumx,&sumy); if (sumx >0) { Zvertex = sumy/sumx+Geo->ver_z_center; ErrorCode = 1; /* vertex found */ } else { ErrorCode = -1; if ( Io->verbose < VERY_LOUD ) cout <<"\t" << "Zvertex : failure (no counts in peak)" << endl; } } /* end if */ else { ErrorCode = -2; if ( Io->verbose < VERY_LOUD ) cout <<"\t" << "Zvertex : failure (peak height =0)" << endl; } if ( ErrorCode <= 0 ) { if ( Io->verbose < VERY_LOUD ) cout << "\t" << " Zvertex : failure ErrorCode =" << ErrorCode; VertexOut[NumOut].ErrorCode = ErrorCode; return STAFCV_OK; /* should really return STAFCV_BAD */ } VertexOut[NumOut].vertex[2] = Zvertex; VertexOut[NumOut].vertexerr[2][2] = ZRESOLUTION*ZRESOLUTION; /* fill diagnostic histogram */ for ( ibin=0; ibinFill(xtmp,z_data[ibin]); } free ( z_data ); /* Z vertex is found, now look for x,y of vertex *************** */ /* basically, we apply same algorithm to find (x,y) vertex position after z position is determined, we assume (x,y) position would be in normal plane which is prependicular to z axis. if we calculate linear equation which is determined from the one pair of inner clump and outer clump we can find a point that meet that normal plane and have the most contents, when we histogram that point . if ( x,y) vertex position is outside the inner barrel, probably ,it is meaningless. from one candidate of outer clump, first find the range that falls into radius of inner barrel. second, find the inner clump that determines linear equation which crosses the normal plane. and then we search for the point that exists in that plane and has the most contents, on the condition that the z vertex position of that point should be in the range ( z_vtx-delz,z_vtx+delz) here delz is determined by how precise z position is resolved. now delz is 0.2cm . in the mean time , we will use two rows of bottom part of barrels to calculate the (x,y) vertex position, which means row number (4 ,5) is used */ x_data=(float *) calloc(Lx,sizeof(float)); if (x_data ==NULL) { if ( Io->verbose < VERY_LOUD ) cout <<"\t" << " X vertex : out of memory ========" << endl; VertexOut[NumOut].ErrorCode=-3; return STAFCV_OK; /* should really be STAFCV_BAD */ } for ( j=0; j 0) { refine_resolution(x_data,Lx,dxstep,xofftmp1,ichmax,&sumx,&sumy); if (sumx > 0 ) { x_vertex=sumy/sumx ; } } /* fill diagnostic histogram */ for ( ibin=0; ibinFill(xtmp,x_data[ibin]); } free(x_data); y_data=(float *) calloc(Ly,sizeof(float)); if (y_data ==NULL) { if ( Io->verbose < VERY_LOUD ) cout <<"\t" << " Y vertex : out of memory ========" << endl; VertexOut[NumOut].ErrorCode=-4; /* no memory for y vertex search */ return STAFCV_OK; /* should really be STAFCV_BAD */ } for ( j=0; j 0) { refine_resolution(y_data,Ly,dystep,yofftmp1,ichmax,&sumx,&sumy); if (sumx > 0 ) { y_vertex=sumy/sumx; } } /* fill diagnostic histogram */ for ( ibin=0; ibinFill(xtmp,y_data[ibin]); } free(y_data); outer_radius=(Geo->voslr+Geo->dvosl[0]); x_vertex=x_vertex-outer_radius; y_vertex=y_vertex-outer_radius; /* using rotation matrix to calculate (x,y ) position in the MVD coordinate system */ Xvertex = (float)cos(30.)*x_vertex+(float)sin(30.)*y_vertex; Yvertex =-(float)sin(30.)*x_vertex+(float)cos(30.)*y_vertex; if ( Io->verbose < LOUD ) cout <<"\t" << "Xvertex= " << Xvertex << endl; if ( Io->verbose < LOUD ) cout <<"\t" << "Yvertex= " << Yvertex << endl; if ( Io->verbose < LOUD ) cout <<"\t" << "Zvertex= " << Zvertex << endl; VertexOut[NumOut].vertex[0] = Xvertex; VertexOut[NumOut].vertex[1] = Yvertex; VertexOut[NumOut].vertexerr[0][0] = dxstep*dxstep; VertexOut[NumOut].vertexerr[1][1] = dystep*dystep; VertexOut[NumOut].ErrorCode=ErrorCode; /*get final cpu time*/ cpu_out=(float)clock()/(float)CLOCKS_PER_SEC; cputime=cpu_out - cpu_in; CpuPerEvP->Fill(cputime); CpuVsOccP->Fill(Cntrl[0].foccup[0],cputime); if ( Io->verbose < LOUD ) cout <<"\t" << "time elapsed=" << cputime << endl; if ( Io->verbose < LOUD ) cout << "}" << endl; return STAFCV_OK; } void find_closest_center(float *data,int Ll,float *max,int *channel) { int ichmax; float ymax; int i; float test1,test2; ichmax=0,ymax=data[0]; for (i=1;i ymax) { ymax=data[i]; ichmax=i; } else { if ( data[i] ==ymax && data[i] >0) { test1=(float)(ichmax-Ll/2)-0.5; test2=(float)(i-Ll/2)-0.5; test1 = (test1>0) ? test1 : -test1; test2 = (test2>0) ? test2 : -test2; if ( test2 < test1) { ymax=data[i]; ichmax=i; } } } } *max=ymax; *channel=ichmax; } void refine_resolution(float *data,int Ll, float step, float z1, int ichmax,float *sumx,float *sumy) /* 30-Jul-98 J.P.Sullivan removed id from calling sequence, */ /* add step, and z1, calculate position corresponding to each */ /* channel from step and z1 rather than HIX call */ { int i,i1,i2; float dzstep; float zval; dzstep=step; i1=ichmax-1; i2=ichmax+1; if (i1 < 0) i1=0; if(i2 > Ll-1) i2=Ll-1; /* Ll-1 is maximum index */ *sumx=0.; *sumy=0.; for (i=i1;i<=i2;++i) { zval= z1 + (((float)(i))+0.5)*step; *sumx=*sumx+data[i]; *sumy=*sumy+data[i]*zval; } } void proj_histogram(DMVDGEO_ST *Geo, TABLE_HEAD_ST *Clmp_h,DMVDCLMP_ST *Clmp , float *data, float x1, float step, int ndim, float z_vtx,int row) /* 30-Jul-98 J.P.Sullivan removed hid from calling sequence and added */ /* data, x1, step, ndim -- the calculations are now based on the data */ /* array rather than the hbook histogram */ { int candidate; float zlower,zupper; int first,last; float dz1,dz2; float proj; int i,j; int nbin; float drtmp,dztmp,drtmp2; int nouter; for ( i=0; i< Clmp_h->nok;++i) if (Clmp[i].shell ==1) { nouter=i; break; } drtmp=(Geo->voslr+Geo->dvosl[0])-(Geo->vislr+Geo->dvisl[0]); drtmp2=(Geo->voslr+Geo->dvosl[0])+(Geo->vislr+Geo->dvisl[0]); for ( i=nouter;i nok;++i) /* if both outer clump and inner clump have the same z vertex position we can not determine the projected position into the normal plane */ if ( Clmp[i].shell== 1 && Clmp[i].row ==row && Clmp[i].zavg != z_vtx) { candidate=i; find_range(Clmp[candidate].zavg,z_vtx,0.2,Geo,&zlower,&zupper); /* calculate inner clump indexes for one outer clump*/ first=last=-1; for ( j=0; j < nouter; ++j) if ( Clmp[j].shell==0 && Clmp[j].row==row){ if (Clmp[j].zavg > zlower && Clmp[j].zavg < zupper &&first == -1) first=j; } for ( j=(nouter-1); j >0; --j) if ( Clmp[j].shell==0 && Clmp[j].row==row){ if (Clmp[j].zavg > zlower &&Clmp[j].zavg < zupper && last ==-1) last=j; } /* if we found inner clumps which falls into the the radius of inner barrel we calculate the position of projected point from the pair of clumps into the normal plane */ dz2=Clmp[candidate].zavg-z_vtx; dz2=dz2 >0 ? dz2 :-dz2; if ( first != -1 ) { for (j=first;j <=last ;++j) { /* impose the condition that the projected z vertex position should be in the range (z_vtx-delz, z_vtx+delz) */ dztmp=(Clmp[j].zavg-Clmp[candidate].zavg)/drtmp*(Geo->voslr+Geo->dvosl[0]) +Clmp[candidate].zavg; if ( dztmp >(z_vtx-0.2) && dztmp <(z_vtx+0.2)) { dz1=Clmp[j].zavg-z_vtx; dz1=dz1>0 ? dz1 : -dz1; proj=((Geo->voslr+Geo->dvosl[0])-(Geo->vislr+Geo->dvisl[0]))/(dz2-dz1)*dz2; /* if ( proj < drtmp || proj > drtmp2 ) cout <<"\t" << "proj : " << proj << endl; */ if ( proj > (drtmp+1) && proj < (drtmp2-1)) { nbin = (int)((proj-x1)/step); if ( (0<=nbin) && (nbinvislr+Geo->dvisl[0]); dr2=(Geo->voslr+Geo->dvosl[0])+(Geo->vislr+Geo->dvisl[0]); if ( z_pos < z_vtx) { dz=z_vtx-z_pos-delz; dztmp=dr1*dz/dr2; *upper=z_vtx+delz; *lower=z_vtx-dztmp-delz; } else { dz=z_pos-z_vtx-delz; dztmp=dr1*dz/dr2; *upper=z_vtx+dztmp+delz; *lower=z_vtx-delz; } }