/*:>-------------------------------------------------------------------- **: FILE: mMvbdNdEta.cc **: HISTORY: **: 00jan93-v000a-hpl- Created by stic Version **: Id: idl.y,v 1.2 1998/02/03 22:21:05 dave Exp **:<------------------------------------------------------------------*/ #include "mMvbdNdEta.h" #include #include #include #include void MvbEta( DMVDDNDETAPAR_ST *dNdEtaPar, DMVDGEO_ST *Geo, DMVBGEO_ST *bGeo, DMVBDBASE_ST *bDbase, DMVDIO_ST *Io, short bin, float z1, float z2, float zvtx, float dphi, float ObsSig, float *EtaAvg, float *dEta, float *EstHit, float *dNdEta, float *dNdEtaErr ); float MvbGeomCorrFactor( DMVDGEO_ST *Geo, float delZ ); #define NO_ERROR 0 /* used for mesage output with no error */ #define VERY_LOUD 4 /* error messages -- hard to turn off */ #define LOUD 3 /* messages which are usually output */ #define QUIET 2 /* used for messages which are output at default setting */ #define SILENT 1 /* for messages which are usually not output */ #define TWOPI 6.283185 /* 2*pi */ #define DELPHI 1.0472 /* 60 degrees (spacing between rows) */ /* translated to radians */ /* these software ID codes are used for the different algorithms used to get */ /* dN/deta/dphi. Whenever a significant change in the algorithm is made, these */ /* should be incremented. */ #define SOFT_ID_WHOLE 1.01 /* software ID code for dN/deta which */ /* combines all rows in one eta bin */ #define SOFT_ID_ROW 2.01 /* software ID code of dN/deta/dphi */ /* which is done for individual rows */ long mMvbdNdEta_( TABLE_HEAD_ST *Geo_h, DMVDGEO_ST *Geo , TABLE_HEAD_ST *bGeo_h, DMVBGEO_ST *bGeo , TABLE_HEAD_ST *Par_h, DMVDPAR_ST *Par , TABLE_HEAD_ST *bPar_h, DMVBPAR_ST *bPar , TABLE_HEAD_ST *bDbase_h, DMVBDBASE_ST *bDbase , TABLE_HEAD_ST *bRaw_h, DMVBRAW_ST *bRaw , TABLE_HEAD_ST *VertexOut_h, DMVDVERTEXOUT_ST *VertexOut , TABLE_HEAD_ST *dNdEtaPar_h, DMVDDNDETAPAR_ST *dNdEtaPar , TABLE_HEAD_ST *dNdEtaOut_h, DMVDDNDETAOUT_ST *dNdEtaOut , TABLE_HEAD_ST *dMultOut_h, DMVDMULTOUT_ST *dMultOut , TABLE_HEAD_ST *Io_h, DMVDIO_ST *Io ) { /*:>-------------------------------------------------------------------- **: ROUTINE : mMvbdNdEta_ **: **: DESCRIPTION : Calculate dN/d(eta) distributions from the MVD using the **: digitized information from the barrel. **: **: AUTHOR : J.H.Park, S.S.Ryu, Y.K.Kim (YonSei Univ.) (Dec 15 97) **: Translated from the original Fortran version ver_neta.f, **: ver_geta.f, ver_gcor.f written by J.P.Sullivan **: **: ARGUMENTS : **: Geo - MVD geometry info which comes from pisa **: Geo_h - header Structure for Geo **: bGeo - contains z of each strip and cell **: bGeo_h - header Structure for bGeo **: bPar - misc constants such as the half length of **: the barrel, the signal per mip in keV,... **: bPar_h - header Structure for bPar **: bRaw - raw data structure for barrel **: bRaw_h - header Structure for bRaw **: VertexOut - vertex finding output table **: VertexOut_h - header structure for VertexOut **: Io - for verbose parameter **: Io_h - header structure for Io **: INOUT: **: dNdEtaPar - misc parameters used to calculate **: dN/deta including the average ADC signal **: per particle (AvgMip) and the number of **: adjacent strips to combine in calculating **: dN/deta (nBunch). **: dNdEtaPar_h - header Structure for dNdEtaPar **: OUT: **: RETURNS: STAF Condition Value **: **: INTERNAL SPECIFICATIONS: **: **: n_ver_zbin = NZDIM/nBunch = number of bins used in dN/d(eta) **: dEta Delta ETA occupied by a group of strips **: dNdEta estimated value of dN/d(eta) from sum **: dNdEta estimated value of dN/d(eta)d(phi) from sum **: EtaAvg = (ETA1+ETA2)/2. **: **: Revision History **: Name Date Comment **: C.F. Maguire Jan. 26, 1998 Initialize ObsSig[48] array to 0 **: **: S.S.Ryu Feb. 10, 1998 changing "atan" in MvbEta to "atan2" **: atan won't give correct value. **: **: J.H.PARK Feb. 10, 1998 delete some arrays(ObsSig,first_ch,..) **: histogram dN/dEta for real vertex **: **: J.P.Sullivan Mar. 5, 1998 set variable vertex when Cntl->choice **: is not either 1 or 2 (Cntl->choice=0 **: means use both algorithms). **: **: J.P.Sullivan Apr. 14, 1998 replace printf with MvdMess calls, **: including in commented out uses of **: printf **: J.P.Sullivan Apr. 28, 1998 remove header table from calling **: sequence. Remove histogram which used **: the true vertex position in calculating **: dN/deta. Add comments and some **: minor code cleanup. **: J.P.Sullivan Jun. 17, 1998 Add dphi correction (the active part of **: each row of Si occupies less than **: 60 degress in phi due to dead areas **: at the edges, etc. Also add dN/deta*dphi. **: Add new staf table dMvddNdEtaOut to **: calling sequence. Define first software **: ID codes: **: SOFT_ID_WHOLE = 1.01 **: SOFT_ID_ROW = 2.01 **: Modified MvbEta calling sequence to **: include dphi and an estimated **: uncertainty on dN/deta. **: Removed histogram filling from this **: module -- it is now done elsewhere **: using information in dMvddNdEtaOut. **: Remove dMvdHst table from calling **: sequence since it is no longer used. **: J.P.Sullivan 10-Aug-98 Remove mMvdPseudoTrk, mMvdZcor from calling **: sequence, add mMvdVertexOut. **: **: S.S.Ryu Sep, 4, 1998 estamate the number of particles making **: use of the calibration data **: **: S.S.Ryu Sep, 8, 1998 Add new staf table dMvdMultOut to **: calling sequence. **: J.P.Sullivan 6-Oct-1998 Add protection for calls with no **: input data (in bRaw). **: J.P.Sullivan 6-Oct-1998 Change to C++ (actually, the code is **: still C). Replace mMvdMess with cout. **: Faster dN/deta algorithm used. **: J.P.Sullivan 12-Oct-1998 Protect against overflows in dNdEtaOut. **: Add code to calculate cpu time used. **: S.S.Ryu 06-Jul-1999 removes unnecessary access to the table **: dMvdTrigControl **: removes unused file stream "FILE *infile" **: changes "static bCall" **: to "static short bCall" ***************************************************************************/ int i,j,k,index; int n_ver_zbin; short NPHDIM=bPar->NPHDIM; short STRIP=Geo->nstrip; short PANEL=bPar->NCHDIM; float vertex; DMVBRAW_ST *btemp; int first_ch; int last_ch; float z_first_ch; float z_last_ch; float dEta,EtaAvg,EstHit,dNdEta,dNdEtadPhi[6]; float dNdEtaErr; /* uncertainty on dN/deta */ float dNdEtadPhiErr[6]; /* uncertainty on dN/deta/dphi */ float ObsSig; float ObsSigPhi[6]; /* Signal for individual phi bins */ float multiplicity=0.; int next; static short bCall=0; float dphi_row; /* dphi (radians) covered by the active part of one row */ /* of the inner barrel. Calculated from geometry table */ float dphi_whole; /* dphi covered by all rows combined 6*DPHI_ROW */ int adc_array_size; // size of adc_sum float *adc_sum; // for array to add up ADC signals int first_error; float cputime,cpu_in,cpu_out; // for tests of cpu time use /****************************************************************************** ******************************************************************************/ cpu_in=(float)clock()/(float)CLOCKS_PER_SEC; bCall++; // count the number of calls first_error = 0; // to prevent numerous similar error messages if ( Io->verbose < LOUD ) cout << "\nmMvbdNdEta{\n\tbCall = " << bCall << endl; /* The following parameters should be set outside this function, however */ /* if they have not been set beforehand, set these default values */ if (dNdEtaPar_h->nok==0){ dNdEtaPar->AvgMip = Par->full_scale/Par->smax_mip; /* ADC channels/mip */ dNdEtaPar->nBunch = 32; /* strips per bunch */ dNdEtaPar->tooClose = 5.0; /* (cm) */ dNdEtaPar_h->nok = 1; } if ( dNdEtaPar->nBunch <= 0 ) { if ( Io->verbose < LOUD ) cout << "\tillegal value of dNdEtaPar->nBunch=" << dNdEtaPar->nBunch << " reset it to 32" << endl; dNdEtaPar->nBunch = 32; // set up a default value to avoid divide by zero } vertex=VertexOut->vertex[2]; // local variable for vertex position /***********************************************/ /* Use the inner barrel to determine dN/d(eta) */ /***********************************************/ n_ver_zbin = bPar->NZDIM/dNdEtaPar->nBunch; /* number of bins in z */ dphi_row = 2.0*(float)atan2( (double)Geo->dvwr1[1], (double)(Geo->vislr+Geo->dvwr1[0]) ); dphi_whole = 6.*dphi_row; if ( Io->verbose < SILENT ) cout << "\n\tn_ver_zbin=" << n_ver_zbin << "\tbRaw_h->nok=" << bRaw_h->nok << "\tdphi_row =" << dphi_row << "\tdphi_whole=" << dphi_whole << endl; // Figure out the size of the "adc" array needed. // It is the number of eta (i.e. z) bins times // the number of phi bins. Then allocate space for // the array. adc_array_size = n_ver_zbin*NPHDIM; adc_sum = (float *)calloc ( adc_array_size, sizeof(float) ); if ( adc_sum == NULL ) { // If the memory could not be allocated, send an error message // and return -- the rest of the module will not work if ( Io->verbose < VERY_LOUD ) cout << "\tmMvbdNdEta out of memory" << endl; return STAFCV_OK; } else { for ( i=0; i < adc_array_size ; i++ ) adc_sum[i]=0; // clear ADC array for ( i=0; i < bRaw_h->nok ; i++ ) { if ( bRaw[i].shell == 0 ) // only the inner shell is used { // the index groups dNdEtaPar->nBunch adjacent channels (in z) index = (bRaw[i].strip + bRaw[i].cell*STRIP + bRaw[i].row*STRIP*PANEL)/ dNdEtaPar->nBunch; if ( (index >= 0) && (index < adc_array_size) ) { adc_sum[index] += bRaw[i].adc; } } } } /* loop is over eta bins (which are based on groups of adjacent strips) */ for ( j=0; jnBunch*j; /* first strip in this bunch */ last_ch = (j+1)*dNdEtaPar->nBunch-1; /* last strip in this bunch */ // should subtract pitch/2 from 1st and add patch/2 to last. // In the first pass I skip this correction since it seems to // have been left out of the previous version and I want to make // sure this version returns the same answer as the previous // version before I fix it. z_first_ch = bGeo->zstrip[first_ch] + Geo->ver_z_center; /* z of 1st strip */ z_last_ch = bGeo->zstrip[last_ch] + Geo->ver_z_center; /* z of last strip */ if ( Io->verbose < SILENT ) cout << "\n\tz_first_ch =" << z_first_ch << "\tz_last_ch =" << z_last_ch << endl; ObsSig = 0.; for (k=0;k< NPHDIM;++k) { index=j + k*n_ver_zbin; if ( index > 0 && index < adc_array_size ) { ObsSig += adc_sum[index]; /* the ADC total for this eta */ ObsSigPhi[k] = adc_sum[index]; /* the ADC total for this phi */ } } if ( Io->verbose < SILENT ) cout << "\tObsSig" << ObsSig << "\n\tObsSigPhi=" << ObsSigPhi[0] << " " << ObsSigPhi[1] << " " << ObsSigPhi[2] << " " << ObsSigPhi[3] << " " << ObsSigPhi[4] << " " << ObsSigPhi[5] << endl; /* calculate dN/deta for all rows combined */ EstHit=0.; MvbEta( dNdEtaPar, Geo, bGeo, bDbase, Io, j, z_first_ch, z_last_ch, vertex, dphi_whole, ObsSig, &EtaAvg, &dEta,&EstHit, &dNdEta, &dNdEtaErr ); if ( dNdEtaOut_h->nok < dNdEtaOut_h->maxlen ) { dNdEtaOut_h->nok++; /* new entry in dN/deta/dphi output table */ next = dNdEtaOut_h->nok; dNdEtaOut[next].id = bCall; /* ID just counts calls to this routine */ dNdEtaOut[next].dndeta = dNdEta; /* dN/deta/dphi value */ dNdEtaOut[next].dndetaerr = dNdEtaErr; /* uncertainty */ dNdEtaOut[next].eta = EtaAvg; /* eta at center of bin */ dNdEtaOut[next].deta = dEta; /* size of eta bin */ dNdEtaOut[next].phi = 0.0; /* phi value = 0 when whole barrel used */ dNdEtaOut[next].dphi = dphi_whole; /* delta phi */ dNdEtaOut[next].softdndeta = SOFT_ID_WHOLE; /* software ID for algorithm using a */ /* all rows at this eta combined */ } else { if ( first_error == 0 ) { if ( Io->verbose < VERY_LOUD ) cout << "\tERROR in mMvbdNdEta -- no more room in dNdEtaOut table" << "\n\tthis entry will be skipped. You should allocate more" << "\n\trows for this table. dNdEtaOut_h->nok=" << dNdEtaOut_h->nok << endl; } first_error = 1; // set to 1 to avoid getting the message above repeatedly } multiplicity+=EstHit; if ( Io->verbose < SILENT ) cout << "\t next=" << next << " EtaAvg=" << EtaAvg << " dEta=" << dEta << " dNdEta=" << dNdEta << " dNdEtaErr" << dNdEtaErr << endl; /* calculate dN/deta for individual rows */ for (k=0;k< NPHDIM;++k) { MvbEta( dNdEtaPar, Geo, bGeo, bDbase, Io, j, z_first_ch, z_last_ch, vertex, dphi_row, ObsSigPhi[k], &EtaAvg, &dEta,&EstHit, &dNdEta, &dNdEtadPhiErr[k] ); dNdEtadPhi[k]=dNdEta; if ( dNdEtaOut_h->nok < dNdEtaOut_h->maxlen ) { dNdEtaOut_h->nok++; /* new entry in dN/deta/dphi output table */ next = dNdEtaOut_h->nok; dNdEtaOut[next].id = bCall; /* ID just counts calls to this routine */ dNdEtaOut[next].dndeta = dNdEtadPhi[k]; /* dN/deta/dphi value */ dNdEtaOut[next].dndetaerr = dNdEtadPhiErr[k]; /* estimated uncertainty on dN/deta */ dNdEtaOut[next].eta = EtaAvg; /* eta at center of bin */ dNdEtaOut[next].deta = dEta; /* size of eta bin */ dNdEtaOut[next].phi = (float)k*DELPHI; /* phi value */ dNdEtaOut[next].dphi = dphi_row; /* delta phi */ dNdEtaOut[next].softdndeta = SOFT_ID_ROW; /* software ID for algorithm using a */ /* single row */ if ( Io->verbose < SILENT ) cout << "\tnext=" << next << " EtaAvg=" << EtaAvg << " dEta=" << dEta << " dNdEtadphi=" << dNdEtadPhi[k] << " dNdEtaErr=" << dNdEtadPhiErr[k] << " phi=" << dNdEtaOut[next].phi << endl; } else { if ( first_error == 0 ) { if ( Io->verbose < VERY_LOUD ) cout << "\tERROR in mMvbdNdEta -- no more room in dNdEtaOut table" << "\n\tthis entry will be skipped. You should allocate more" << "\n\trows for this table. dNdEtaOut_h->nok=" << dNdEtaOut_h->nok << endl; } first_error=1; // set to 1 to avoid getting the message above repeatedly } } } dMultOut->id =bCall; dMultOut->totalmult =multiplicity; dMultOut->totalmulterr=(float)sqrt((double)multiplicity); dMultOut->softmult =SOFT_ID_WHOLE; dMultOut_h->nok=1; free ( adc_sum ); /*get final cpu time*/ cpu_out=(float)clock()/(float)CLOCKS_PER_SEC; cputime=cpu_out - cpu_in; if ( Io->verbose < LOUD ) cout << "\tcputime=" << cputime << endl; if ( Io->verbose < LOUD ) cout << "}" << endl; return STAFCV_OK; } void MvbEta( DMVDDNDETAPAR_ST *dNdEtaPar, DMVDGEO_ST *Geo, DMVBGEO_ST *bGeo, DMVBDBASE_ST *bDbase, DMVDIO_ST *Io, short bin, float z1, float z2, float zvtx, float dphi, float ObsSig, float *EtaAvg, float *dEta, float *EstHit, float *dNdEta, float *dNdEtaErr ){ /*: DESCRIPTION: calculate dN/deta for a single bin **: **: ARGUEMENTS to MvbEta: **: **: IN: dNdEtaPar -- misc parameters used to calculate **: dN/deta including the average ADC signal **: per particle (AvgMip) and the number of **: adjacent strips to combine in calculating **: dN/deta (nBunch). A staf table. **: Geo -- MVD geometry info which comes from pisa. **: A staf table. **: z1 -- z (cm) of first strip in this bunch **: z2 -- z (cm) of last strip in this bunch **: zvtx -- z of vertex (cm) **: dphi -- delta phi (azimuthal angle) coverage of **: this bin (radians) **: ObsSig -- total ADC signal from all strips in the bunch **: Io -- for verbose parameters **: **: OUT: *EtaAvg -- average eta in bin =(eta(z1)+eta(z2))/2 **: *dEta -- deta of this bin = eta(z2)-eta(z1) **: *EstHit -- estimated number of charged particles which **: hit this bunch of strips **: *dNdEta -- dN/deta, where dN is calculated from ObsSig, **: misc geometry parameters, and the average **: signal per mip. **: *dNdEtaErr-- uncertainty on dN/deta ************************************************************************/ long index; float MiddleZ; /* Z at middle of bunch (cm) */ float delZ; /* MiddleZ - (z of vertex) (cm) */ float CorrFactor; /* geometry correction factor -- 1 for a */ /* particle incident normal to the surface */ /* of the Si panel. */ float AvgSig; /* average signal per particle (ADC channels) */ /* for a particle with the expected angle of */ /* incidence */ float theta1; /* theta (polar angle) corresponding to z1 */ float theta2; /* theta (polar angle) corresponding to z2 */ float tanTmp; /* intermediate variable in eta calculation */ float Eta1; /* eta corresponding to z1 */ float Eta2; /* eta corresponding to z2 */ double a,b,c,d,f; /* ********************************************************************* */ MiddleZ = (z1+z2)/2.; /* middle of bunch */ delZ = MiddleZ-zvtx; /* dz of bunch */ if ( Io->verbose < SILENT ) cout << "\tz1=" << z1 << " z2=" << z2 << " MiddleZ=" << MiddleZ << " delZ=" << delZ << endl; if ( Io->verbose < SILENT ) cout << "\tCorrFactor=" << CorrFactor << " AvgSig= " << AvgSig << " EstHit =" << EstHit << endl; theta1 = (float)atan2( (double)Geo->vislr,(double)(z1-zvtx) ); theta2 = (float)atan2( (double)Geo->vislr,(double)(z2-zvtx) ); tanTmp = (float)tan( (double)(theta1/2.)); Eta1 = -(float)log( (double)tanTmp ); /* eta(z1) */ if ( Io->verbose < SILENT ) cout << "\tGeo->vislr=" << Geo->vislr << " z1-zvtx =" << z1-zvtx << " theta1 = " << theta1 << " theta2 = " << theta2 << " tanTmp = " << tanTmp << " Eta1 = " << Eta1 << endl; tanTmp = (float)tan( (double)(theta2/2.)); Eta2 = -(float)log( (double)tanTmp ); /*eta(z2) */ if ( Io->verbose < SILENT ) cout << "\ttanTmp=" << tanTmp << " Eta2=" << Eta2 << endl; /* ******************************************** */ /* now setup the variables which are returned */ *EtaAvg = (Eta1+Eta2)/2.; *dEta = Eta2-Eta1; index=dNdEtaPar->nBunch*bin; a=bDbase[index].calib_const[0]; b=bDbase[index].calib_const[1]; c=bDbase[index].calib_const[2]; d=bDbase[index].calib_const[3]; f=bDbase[index].calib_const[4]; AvgSig = a + b*(*EtaAvg)*(*EtaAvg) + c*(*EtaAvg)*(*EtaAvg)*(*EtaAvg)*(*EtaAvg) + exp( d + f*fabs((double)*EtaAvg) ); CorrFactor = MvbGeomCorrFactor( Geo, delZ); /* geometry correction */ *EstHit=ObsSig/(AvgSig*CorrFactor); *dNdEta =*EstHit/((*dEta)*(dphi/TWOPI)); if ( *EstHit>0. ) { *dNdEtaErr= *dNdEta/(float)sqrt((double)(*EstHit)); } else { *dNdEtaErr=0.; } if ( Io->verbose < SILENT ) cout << "\tEtaAvg= " << *EtaAvg << " dEta= " << *dEta << " dNdEta= " << *dNdEta << " dphi= " << dphi << endl; } float MvbGeomCorrFactor( DMVDGEO_ST *Geo, float delZ ){ /* **: DESCRIPTION: Calculate geometry correction factor. The calculation **: is of 1/cos(alpha), where alpha is the angle of incidence **: of a particle hitting the MVD inner barrel at a distance **: delZ from the vertex. The angle is defined such that alpha=0 **: is a particle entering normal to the Si. If the finite **: azimuthal coverage of each strip was ignored, then this **: angle could be trivially calculated (it would be sin(theta), **: where theta is the polar angle of the particle). However, **: this calculation DOES consider the finite azimuthal coverage **: of each Si panel in the barrel. An average over a uniform **: distribution in azimuthal angles is folded into the polar **: angle to get an average angle of incidence. The integral is **: too messy to try to type into a comment in a program. It is **: shown in an entry in my (JPSullivan) logbook dated 12-Jan-93. **: **: ARGUEMENTS to MvbEta: **: **: IN: Geo -- MVD geometry info which comes from pisa. **: A staf table. **: delz -- (z of center of bunch) - (z of vertex) (cm) **: **: RETURNS: Geometry correction factor (1/cos(alpha) -- details above) **: ************************************************************************/ float R; /* R of barrel (it has a hexagonal cross section */ /* so this is the closest distance from the */ /* to the beamline) (cm) */ float wid; /* half width of a Si pannel (azimuthal) (cm) */ float sinMax; /* sin of angle (azimuthal direction) between a */ /* vector from the beamline to the center of a Si*/ /* panel and a vector from the beamline to the */ /* end of the same panel */ float phiMax; /* angle corresponding to sinMax */ float absZ; /* absolute value of delz */ float tan_theta; /* tan polar angle to center of bunch of strips */ float cos_theta; /* cos polar angle to center of bunch of strips */ /* ******************************************************************** */ R = Geo->vislr; /* local copy of the barrel radius */ wid = Geo->dvisl[1]; /* local copy of half-width of panel */ sinMax = wid/(float)sqrt( (double)(R*R+wid*wid) ); phiMax = (float)asin( (double)sinMax ); absZ = (float)fabs((double)delZ); /* These is a special case for small absZ where there is a potential */ /* 0/0 in the calculation, for that special case use this simpler */ /* limit of the general formula. Otherwise do the normal ... */ if ( absZ <= 0.1) return phiMax/sinMax; else { tan_theta= R/absZ; /* cos, then tan of polar angle theta */ cos_theta = absZ/(float)sqrt( (double)(absZ*absZ+R*R) ); return phiMax/( tan_theta*(float)asin( (double)(cos_theta*sinMax) ) ); } }