/*:>-------------------------------------------------------------------- **: FILE: mMvcdNdEta.cc **: HISTORY: **: 00jan93-v000a-hpl- Created by stic Version **: Id: idl.y,v 1.2 1998/02/03 22:21:05 dave Exp **:<------------------------------------------------------------------*/ #include "mMvcdNdEta.h" #include #include "PhHistogramFactory.hh" #include float MvcNum( int iADCCor, float xnBar, int nMax, float hn[][256], float gn[][256] ); void MvcEta( DMVDDNDETAPAR_ST *dNdEtaPar, float r1, float r2, float zvtx, int nh_Ring, float z, float *dEta, float *EtaAvg, float *dNdEta, float *dNdEtaErr ); #define MaxIter 2 /* max iterations */ #define ETA_MIN_HIST -7. /* min eta in histograms */ #define ETA_MAX_HIST 7. /* max eta in histograms */ #define N_ETA_BINS_IN_HIST 56 /* number of eta bins in histograms */ #define TWOPI 6.283185 #define SOFT_ID_AVGDEDX 3.01 /* software ID code for dN/deta from */ /* average dE/dx per particle */ #define SOFT_ID_DECONV 4.01 /* software ID code for dN/deta from */ /* deconvolution method */ #define LOUD 3 #define VERY_LOUD 4 // intended for error messages long mMvcdNdEta_( TABLE_HEAD_ST *Geo_h, DMVDGEO_ST *Geo , TABLE_HEAD_ST *cGeo_h, DMVCGEO_ST *cGeo , TABLE_HEAD_ST *Par_h, DMVDPAR_ST *Par , TABLE_HEAD_ST *cPar_h, DMVCPAR_ST *cPar , TABLE_HEAD_ST *cRaw_h, DMVCRAW_ST *cRaw , TABLE_HEAD_ST *VertexOut_h, DMVDVERTEXOUT_ST *VertexOut , TABLE_HEAD_ST *Cntrl_h,DMVDTRIGCONTROL_ST *Cntrl , TABLE_HEAD_ST *dNdEtaPar_h, DMVDDNDETAPAR_ST *dNdEtaPar , TABLE_HEAD_ST *dNdn_h, DMVDDNDETA_ST *dNdn , TABLE_HEAD_ST *dNdEtaOut_h, DMVDDNDETAOUT_ST *dNdEtaOut , 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 pads. **: Distributions which use the real vertex and **: which uses the vertex found by the tracking algorithm **: are filled. This gets called once per event. **: **: 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_peta.f,vpad_num.f written by J.P.Sullivan **: **: ARGUMENTS: **: IN: **: Geo - MVD geometry parameters from pisa **: Geo_h - header Structure for Geo **: cGeo - geometry constants related to pads **: cGeo_h - header Structure for cGeo **: Par - misc parameters such as energy loss **: per mip, ... **: Par_h - header Structure for Par **: cPar - misc dimensions of pad detectors **: cPar_h - header Structure for cPar **: cRaw - raw data structure for pads **: cRaw_h - header Structure for cRaw **: VertexOut - vertex finding output table **: VertexOut_h - header Structure for VertexOut **: Cntrl - control parameters, including the **: the parameter (choice) which **: tell which vertex finding algorithm **: was used for this event. **: Cntrl_h - header Structure for Cntrl **: Io - contains verbose parameter **: Io_h - header structure for Io **: INOUT: **: dNdEtaPar - some key parameters for the dN/deta **: calculation including AvgMip and **: tooClose **: dNdEtaPar_h - header Structure for dNdEtaPar **: OUT: **: dNdn - number of hits per pad detector **: dNdn_h - header Structure for dNdn **: **: RETURNS: STAF Condition Value **: **: AvgMip = AVeraGe signal of a Minimum Ionizing Particle at normal **: incidence angle in chan. This should be about 32chan, but **: it might be necessary to adjust this parameter upwards **: a little since the average energy loss should be greater **: than the average for a mip. The value is calculated from **: the scale parameters given in VERPARS. **: tooClose is a parameter which prevents the calculation of **: dN/deta from a pad plane when the calculated vertex **: is too close to the pad plane. **: **: INTERNAL SPECIFICATIONS: **: **: ObsSig the OBServed SIGnal (chan) in a given group of pads. **: dEta delta ETA occupied by pads on a ring **: dNdEta estimated value of dN/d(eta) from sum **: EtaAvg = (ETA1+ETA2)/2. **: **: variables asociated with the pad deconvolution code **: **: fn[n][i] = probability that an event with n+1 particles **: hitting a pad will give and ADC value of i **: gn[n][i] = fn[n][i]/(n+1)! **: hn[n][i] = gn*(n+1) **: **: Revision History **: Name Date Comment **: J.H.Park Feb. 10, 1998 delete some arrays(ObsSig,nh_Ring) **: Histogram dN/dEta for real vertex **: change "tan" to "atan2" **: **: S.S.Ryu 13-Feb.-1998 fixing gn[n][i] and hn[n][i]. **: **: J.P.Sullivan Mar. 6, 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 calls with MvdMess_ **: calls, including those in comments. **: Remove dMvcGeo from calling sequence **: to MvcEta (since it was not used). **: **: J.P.Sullivan Apr. 29, 1998 Remove header from calling sequence, **: delete code which used true vertex **: position. Add/change some comments. **: Change nh4 to nh1 and nh6 to nh2. **: Remove nh3, nh5 which were previously **: the histograms which used the true **: vertex position. Modify MvcEta. **: Fixed code to read in ver_calib.dat. **: J.P.Sullivan June. 17, 1998 Add code to save dN/deta in an output **: structure (dMvddNdEtaOut). Changed **: calling sequence to MvcEta to add a **: calculation of the uncertainty. **: J.P.Sullivan 4-Aug-98 Remove dMvdPseudoTrk and dMvdZcor tables from **: calling sequence, add mMvdVertexOut. Change **: file name to .cc to force compilation with C++ **: compiler. Replace mMvdMess with cout. Replace **: hbook with phenix generic histogram interface. **: Remove dMvdHst from calling sequence. **: J.P.Sullivan 23-Sep-98 Add dMvdIo to calling sequence and use it to **: control volume of output **: **: J.P.Sullivan 14-Oct-98 Added protection to avoid overflowing the **: dNdEtaOut output table. Add code to calculate **: CPU time used. Change algorithm to make it **: faster. **: S.S.Ryu 06-Jul-99 changes the variable "Cntl" to "Cntrl" to make **: this variable consistent throughout the MVD **: codes. ***************************************************************************/ int nMax=4; int MvdMaxADC=Par->maxadc_ver; int ih_test[13]; float xmin,xmax; char ver_calib[]="ver_calib.dat"; FILE *in; float adc[4][256]; static float fn[4][256]; static float gn[4][256]; static float hn[4][256]; float sum1; float xNCalc,xnBar; int factorial; int row,end; float ObsSig; float ErrBar[256]; float dN; int nh_Ring; int IADCCor; float CosRing; float PadDeconv[2]; float xtemp; float dEta,EtaAvg,dNdEta,dNdEta2; float dNdEtaErr,dNdEtaErr2; int next; int iter; int first_error; // flag to avoid repeating err. mess. float fdum; static short cCall=0; int i,j,k,n; float vertex; PhHistogramFactory *hf=PhHistogramFactory::instance(); static PhH1 *Input1Part; static PhH1 *Input2Part; static PhH1 *Input3Part; static PhH1 *Input4Part; static PhH1 *Calc1Part; static PhH1 *Calc2Part; static PhH1 *Calc3Part; static PhH1 *Calc4Part; static PhH1 *NcalcVsADCocc1; static PhH1 *NcalcVsADCocc10; static PhH1 *NcalcVsADCocc25; static PhH1 *NcalcVsADCocc50; static PhH1 *NcalcVsADCocc100; /* the following ADC distribution comes from a staf */ /* calculation using pisa events as input. It is used as*/ /* the ADC distribution expected for a single particle */ /* hitting a single pad. Actually, it is for events with*/ /* average occupancy in the pads less than 0.01, which */ /* means the double hit probability is low, but not */ /* zero. The code atempts to get this distribution from */ /* a data file, but uses this distribution as a default */ /* if the file is missing. */ int iadc1[256] ={ 81, 57, 27, 12, 13, 23, 62, 381, 332, 402, 335, 577, 661, 977, 1461, 2201, 3270, 4028, 4134, 3864, 3045, 2459, 1821, 1396, 1125, 937, 844, 635, 549, 483, 447, 394, 322, 342, 260, 245, 232, 208, 200, 190, 191, 150, 169, 138, 132, 142, 131, 117, 96, 101, 90, 86, 68, 75, 82, 64, 58, 57, 72, 54, 55, 46, 53, 44, 40, 39, 39, 36, 39, 41, 39, 35, 37, 42, 45, 32, 27, 27, 29, 27, 25, 27, 16, 24, 13, 25, 28, 17, 21, 17, 19, 16, 14, 19, 25, 15, 21, 22, 13, 17, 13, 19, 18, 14, 16, 16, 6, 10, 14, 13, 10, 16, 9, 8, 8, 11, 13, 8, 8, 7, 7, 8, 5, 16, 6, 11, 5, 11, 9, 7, 3, 9, 15, 5, 6, 10, 5, 3, 5, 3, 9, 5, 3, 7, 6, 7, 3, 8, 4, 5, 1, 8, 6, 3, 6, 8, 5, 4, 4, 9, 6, 7, 5, 5, 1, 5, 3, 6, 8, 4, 8, 3, 2, 5, 7, 2, 4, 3, 9, 4, 5, 9, 2, 6, 4, 2, 3, 1, 5, 6, 5, 2, 6, 3, 5, 1, 2, 2, 6, 2, 3, 2, 2, 1, 5, 2, 2, 4, 2, 2, 2, 2, 1, 3, 2, 2, 7, 4, 3, 6, 5, 1, 2, 5, 2, 4, 2, 2, 4, 4, 2, 2, 0, 3, 5, 2, 4, 4, 4, 1, 7, 4, 7, 8, 11, 7, 9, 10, 14, 24, 22, 29, 43, 41, 8, 0}; /* ****************************************************** */ /* Calculate a dN/d(eta) distribution based on the energy */ /* deposited in each pad. */ if ( Io->verbose < LOUD ) cout << "\nmMvcdNdEta{" << endl; first_error = 0; /* This table should have been filled before calling this */ /* function, but if it is not, put in some default values */ if ( dNdEtaPar_h->nok == 0) { dNdEtaPar->AvgMip = Par->full_scale/Par->smax_mip; dNdEtaPar->nBunch = 32; /* not used in this routine */ /* but initialize it when */ /* initializing this table */ dNdEtaPar->tooClose = 5.0; /* cm */ dNdEtaPar_h->nok=1; } if ( cCall==0 ) { /* initialize code for deconvolution of pad information */ /* diagnostic histograms for the pad unfolding routines */ xmin=0.; xmax=(float)(MvdMaxADC); Input1Part = hf->CreateH1F ( "Input1Part", "input 1 part dist", MvdMaxADC,xmin,xmax); Input2Part = hf->CreateH1F ( "Input2Part", "input 2 part dist", MvdMaxADC,xmin,xmax); Input3Part = hf->CreateH1F ( "Input3Part", "input 3 part dist", MvdMaxADC,xmin,xmax); Input4Part = hf->CreateH1F ( "Input4Part", "input 4 part dist", MvdMaxADC,xmin,xmax); Calc1Part = hf->CreateH1F ( "Calc1Part", "calc. 1 part dist (bootstrap)", MvdMaxADC,xmin,xmax); Calc2Part = hf->CreateH1F ( "Calc2Part", "calc. 2 part dist (bootstrap)", MvdMaxADC,xmin,xmax); Calc3Part = hf->CreateH1F ( "Calc3Part", "calc. 3 part dist (bootstrap)", MvdMaxADC,xmin,xmax); Calc4Part = hf->CreateH1F ( "Calc4Part", "calc. 4 part dist (bootstrap)", MvdMaxADC,xmin,xmax); NcalcVsADCocc1 = hf->CreateH1F ( "NcalcVsADCocc1", "Ncalc vs ADC, 0.01", MvdMaxADC,xmin,xmax); NcalcVsADCocc10 = hf->CreateH1F ( "NcalcVsADCocc10", "Ncalc vs ADC, 0.10", MvdMaxADC,xmin,xmax); NcalcVsADCocc25 = hf->CreateH1F ( "NcalcVsADCocc25", "Ncalc vs ADC, 0.25", MvdMaxADC,xmin,xmax); NcalcVsADCocc50 = hf->CreateH1F ( "NcalcVsADCocc50", "Ncalc vs ADC, 0.50", MvdMaxADC,xmin,xmax); NcalcVsADCocc100 = hf->CreateH1F ( "NcalcVsADCocc100","Ncalc vs ADC, 1.00", MvdMaxADC,xmin,xmax); in= fopen ("ver_calib.dat","r"); if( in !=NULL ) { if ( Io->verbose < LOUD ) cout << "\tver_calib.dat opened" << endl; while(!feof(in)) { for (j=0; j<=MvdMaxADC; ++j) { fscanf(in,"%f %f %f %f %f",&fdum,&adc[0][j],&adc[1][j],&adc[2][j],&adc[3][j]); } } if ( Io->verbose < LOUD ) cout << "\tver_calib.dat read" << endl; } else { /* only the 1 particle distribution is set here -- there other */ /* distributions are in the data file but not stored locally */ /* in this subroutine. The 2,3,4 particle distributions in */ /* this array are actually only used to check whether the */ /* folded 1 particle distribution, which is used to construct */ /* distributions for two or more hits, agrees with them. Here, */ /* to avoid many large local data arrays this is not done. */ for(i=0; i<=MvdMaxADC; ++i) { adc[0][i]=(float)iadc1[i]; adc[1][i]=0.; adc[2][i]=0.; adc[3][i]=0.; } } /* work with normalized values (probabilities) */ for ( n=0; nFill ((float)i,fn[n][i]); else if ( n==1 ) Input2Part->Fill ((float)i,fn[n][i]); else if ( n==2 ) Input3Part->Fill ((float)i,fn[n][i]); else if ( n==3 ) Input4Part->Fill ((float)i,fn[n][i]); } } } for ( i=0; i<=MvdMaxADC; i++) Calc1Part->Fill((float)i, fn[0][i]); /* estimate 2 ~ nMax particle probability based on the 1 particle */ /* probability distribution. This is a sort of booststrap method. */ for (n=1; nFill((float)i, fn[n][i]); else if ( n==2 ) Calc3Part->Fill((float)i, fn[n][i]); else if ( n==3 ) Calc4Part->Fill((float)i, fn[n][i]); } /* calculate the quantities which are more useful in estimating */ /* the number of hits */ factorial=1; for (n=0; nFill( (float)i,xNCalc); else if ( j == 1 ) NcalcVsADCocc10->Fill( (float)i,xNCalc); else if ( j == 2 ) NcalcVsADCocc25->Fill( (float)i,xNCalc); else if ( j == 3 ) NcalcVsADCocc50->Fill( (float)i,xNCalc); else if ( j == 4 ) NcalcVsADCocc100->Fill((float)i,xNCalc); } } if ( Io->verbose < LOUD ) cout << "\tend of first call section" << endl; } /* for cCall==0 */ cCall++; /* count call to this routine */ /* now we get to the part which happens for every event */ /* when the calculated vertex is close to the pads plane, the dN/dEta */ /* calculation will not work very well--so skip the calculation in */ /* that case. */ /* make a local copy of the vertex position found. The algorithm used */ /* is controlled by Cntrl->choice. */ if ( Cntrl->choice ==1) { vertex=VertexOut->vertex[2]; } else if (Cntrl->choice ==2) { vertex=VertexOut->vertex[2]; } else { /* choose the correlation method Cntrl->choice is neither 1 nor 2 */ vertex=VertexOut->vertex[2]; } for ( k=0; k<2; k++) { /* loop over both pad planes */ // cout << "\tvertex=" << vertex << " cGeo->pad_z[k]=" << cGeo->pad_z[k] // << " dNdEtaPar->tooClose=" << dNdEtaPar->tooClose << endl; if ( (float)fabs((double)(vertex-cGeo->pad_z[k])) >= dNdEtaPar->tooClose ){ // cout << "\tnot too close" << endl; PadDeconv[k]=0; /* add up total signal and count the number of pads "on" */ /* loop over radial bins(=eta bins) */ for ( j=0; jNVPAD_RBINS; ++j) { ObsSig =0.; nh_Ring=0; for ( i=0; inok; ++i) if ( cRaw[i].end==k && cRaw[i].row ==j && cRaw[i].adc>0 ){ ObsSig+=cRaw[i].adc; nh_Ring++; } // cout << "\tj(ring)=" << j << " k(end)= " << k // << " rcuts=" << cGeo->rcuts[j] << " , " << cGeo->rcuts[j+1] // << " z of pad= " << cGeo->pad_z[k] << " ObsSig=" // << ObsSig << " nh_Ring=" << nh_Ring << " calc vertex=" // << vertex << endl; MvcEta(dNdEtaPar, cGeo->rcuts[j], cGeo->rcuts[j+1], vertex, nh_Ring, cGeo->pad_z[k], &dEta, &EtaAvg, &dNdEta, &dNdEtaErr ); if ( dNdEtaOut_h->nok < dNdEtaOut_h->maxlen ) { dNdEtaOut_h->nok++; /* new entry in dN/deta output table */ next = dNdEtaOut_h->nok; dNdEtaOut[next].id = cCall; /* 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 = TWOPI; /* delta phi */ dNdEtaOut[next].softdndeta = SOFT_ID_AVGDEDX; /* software ID for algorithm using avg */ /* dE/dx */ // cout << "\tpad count method next=" << next // << " EtaAvg= " << EtaAvg << " dNdEta=" << dNdEta << 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 } if (nh_Ring>0){ xtemp=cGeo->rvals[j]/(vertex-cGeo->pad_z[k]); CosRing=(float)1./sqrt((double)(1.+xtemp*xtemp) ); for ( iter=0; iterNVPAD_PHI; else xnBar=(float)dN/(float)cPar->NVPAD_PHI; dN=0; for ( i=0; inok; ++i) if ( cRaw[i].row ==j && cRaw[i].end==k ){ IADCCor=(short)ceil((double)cRaw[i].adc*CosRing); xNCalc=MvcNum( IADCCor, xnBar, nMax, hn, gn ); dN+=xNCalc; } } /* for iter */ PadDeconv[k]+=dN; dNdEta2=dN/dEta; if ( dN > 0.0 ) dNdEtaErr2 = dNdEta2/(float)sqrt((double)dN); else dNdEtaErr2 = 0; if ( dNdEtaOut_h->nok < dNdEtaOut_h->maxlen ) { dNdEtaOut_h->nok++; /* new entry in dN/deta output table */ next = dNdEtaOut_h->nok; dNdEtaOut[next].id = cCall; /* ID just counts calls to this routine */ dNdEtaOut[next].dndeta = dNdEta2; /* dN/deta value */ dNdEtaOut[next].dndetaerr = dNdEtaErr2; /* 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 = TWOPI; /* delta phi */ dNdEtaOut[next].softdndeta = SOFT_ID_DECONV; /* software ID for deconvolution alg. */ // cout << "\tPoisson method next=" << next << " EtaAvg=" << EtaAvg // << " dNdEta=" << dNdEta2 << " dEta=" << dEta << 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 } } else{ dN=0; dNdEta2=0; dNdEtaErr2=0.; // cout << "\tcalc vtx, Poisson EtaAvg=" << EtaAvg // << " dNdEta=" << dNdEta2 << endl; if ( dNdEtaOut_h->nok < dNdEtaOut_h->maxlen ) { dNdEtaOut_h->nok++; /* new entry in dN/deta output table */ next = dNdEtaOut_h->nok; dNdEtaOut[next].id = cCall; /* ID just counts calls to this routine */ dNdEtaOut[next].dndeta = dNdEta2; /* dN/deta value */ dNdEtaOut[next].dndetaerr = dNdEtaErr2; /* 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 = TWOPI; /* delta phi */ dNdEtaOut[next].softdndeta = SOFT_ID_DECONV; /* software ID for deconvolution alg. */ } 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 } } } /* for j(radial bins) */ } } dNdn->PadDeconv[0] = PadDeconv[0]; dNdn->PadDeconv[1] = PadDeconv[1]; dNdn_h->nok=1; if ( Io->verbose < LOUD ) cout << "}" << endl; return STAFCV_OK; } float MvcNum( int iADCCor, float xnBar, int nMax, float hn[][256], float gn[][256] ){ int n; float sum1=0.,sum2=0.; float xnBarN=1.; for (n=0; n0) return sum1/sum2; else return 0.; } void MvcEta( DMVDDNDETAPAR_ST *dNdEtaPar, float r1, float r2, float zvtx, int nh_Ring, float z, float *dEta, float *EtaAvg, float *dNdEta, float *dNdEtaErr ){ /*:>----------------------------------------------------------------- **: ROUTINE: MvcEta (in file mMvcdNdEta.c) **: DESCRIPTION: Calculate size of bin in pseudorapidity (*dEta), **: center of bin in pseudorapidity (*EtaAvg), dN/deta **: (*dN/Eta) from the input geometry variables for **: one ring of the pad detector (r1, r2, zvtx, z), **: the number of pads hit in one ring of the pad **: detector (nh_Ring). **: AUTHOR: J.J.Park, S.S.Ryu, Y.K.Kim (YonSei Univ.) (Dec 15 97) **: Translated from the original fortran version which **: was written by J.P.Sullivan. **: ARGUEMENTS: **: IN: *dNdEtaPar misc parameters for dN/deta calc. **: r1 inner radius of this ring of pads (cm) **: r2 outer radius of this ring of pads (cm) **: zvtx z vertex position (cm) **: nh_Ring number of pads hit in this ring **: OUT: *dEta size of pseudorapidity bin formed by **: this ring of pads **: *EtaAvg center of pseudorapidity bin **: *dNdEta dN/deta **: RETURNS: (void) **: INTERNAL SPECIFICATIONS: **: **: REVISION HISTORY: **: Name Date Comment **: J.P.Sullivan 14-Apr-98 deleted dMvcGeo table from input **: since it was not used **: J.P.Sullivan 29-Apr-98 I realized that the dN/deta **: algorithm I was using did not work **: when the crosstalk is non-zero. New **: algorithm is based on number of hits **: in this ring, not the ADC sum. **: nh_Ring replaces ObsSig in calling **: sequence. ***************************************************************************/ float MiddleR,delZ; float theta1,theta2, tanTmp, Eta1,Eta2; MiddleR = (r1+r2)/2.; delZ = z-zvtx; // cout << "\tdelZ=" << delZ << " z= " << z << " zvtx" << zvtx << endl; // cout << "\tnh_Ring = " << nh_Ring << endl; theta1 = (float)atan2( (double)r1,(double)delZ ); tanTmp = (float)tan( (double)theta1/2. ); Eta1 = -(float)log( (double)tanTmp ); theta2 = (float)atan2( (double)r2,(double)delZ ); tanTmp = (float)tan( (double)theta2/2. ); Eta2 = -(float)log( (double)tanTmp ); // cout << "\tEta1=" << Eta1 << " \tEta2=" << Eta2 << endl; *EtaAvg = (Eta1+Eta2)/2.; *dEta = (float)fabs( (double)(Eta2-Eta1) ); *dNdEta = (float)nh_Ring/(*dEta); if ( nh_Ring > 0 ) { *dNdEtaErr= *dNdEta/(float)sqrt((double)nh_Ring); } else { *dNdEtaErr=0.; } // cout << "\tEtaAvg=" << *EtaAvg << " \tdEta=" << *dEta // << " \tdNdEta=" << *dNdEta << endl; return; }