if ( ztest >= maxh && ( ( m_applyWeightThreshold && ztest>m_weightThreshold ) || !m_applyWeightThreshold ) ) { apply threshold later - should we wish it
if we are in triplet mode, even a single pair means 3 consistent hits also bomb out if no maximum (above threshold) is found
at this point we have the histogram with the highest N vertices removed so we can find the "non vertex" pedestal from these, although it will be somewhat lower than perhaps it should be, in case some of the "vertices" we are removing are just random upwards fluctuations
NB: have moved pedestal calculation to before the extraction of the vertices if we calculate it after, then we have too low a pedestal if some vertices are really random fluctuations. If we calculate it before then we overestimate the pedestal, really we should try to decide how many real vertices we have, and then only exclude them, but that level of detail is probably not justified by the correlation with the offline track multiplicity on the vertex
copy vertices to output vector - this is so we can first impose cuts on the vertices we have found should we wish to
528 std::vector<vertex>*
output =
new std::vector<vertex>();
530 long nsp = spVec.size();
531 if ( !nsp )
return output;
539 std::vector<double>
phi(nsp),
rho(nsp), zed(nsp);
555 const double invZBinSize = 1./ZBinSize;
558 const long HalfZhistoBins = long ( ceil( 0.5*(zMax-zMin)*invZBinSize ) );
559 const long NumZhistoBins = 2*HalfZhistoBins;
580 std::vector< std::vector<double> >* zHisto =
m_zHisto;
589 for (
unsigned int sliceIndex = 0; sliceIndex <
m_NumPhiSlices; sliceIndex++ )
596 int allSlicesSize = allSlices.size();
598 for (
long pointIndex = 0; pointIndex < nsp; pointIndex++ )
604 allSlices[
phiIndex ]->AddPoint(
rho[ pointIndex ], zed[ pointIndex ],
phi[ pointIndex ], lyr[ pointIndex ] );
608 std::vector< std::vector< double > > allLayerRhos, allLayerZs, allLayerPhis;
613 for (
unsigned int sliceIndex = 0; sliceIndex <
m_NumPhiSlices; sliceIndex++ )
615 allSlices[ sliceIndex ]->MakeWideLayers( &allLayerRhos, &allLayerZs, &allLayerPhis, &allSliceWidths,
filled, &filledLayers );
622 for (
unsigned int sliceIndex = 0; sliceIndex <
m_NumPhiSlices; sliceIndex++ )
624 nHisto[0][sliceIndex].resize( NumZhistoBins + 1 );
625 zHisto[0][sliceIndex].resize( NumZhistoBins + 1 );
629 for (
unsigned int sliceIndex = 0; sliceIndex <
m_NumPhiSlices; sliceIndex++ )
631 nHisto[1][sliceIndex].resize( NumZhistoBins + 1 );
632 zHisto[1][sliceIndex].resize( NumZhistoBins + 1 );
637 for (
unsigned int sliceIndex = 0; sliceIndex <
m_NumPhiSlices; sliceIndex++ )
639 allSlices[ sliceIndex ]->GetHistogram( &(
nHisto[0][sliceIndex] ), &( zHisto[0][sliceIndex] ),
640 &(
nHisto[1][sliceIndex] ), &( zHisto[1][sliceIndex] ),
645 delete allSlices[ sliceIndex ];
651 nHisto[0][0].resize( NumZhistoBins + 1 );
652 zHisto[0][0].resize( NumZhistoBins + 1 );
655 nHisto[1][0].resize( NumZhistoBins + 1 );
656 zHisto[1][0].resize( NumZhistoBins + 1 );
660 for (
unsigned int sliceIndex = 0; sliceIndex <
m_NumPhiSlices; sliceIndex++ )
662 allSlices[ sliceIndex ]->GetHistogram( &(
nHisto[0][0] ), &( zHisto[0][0] ),
663 &(
nHisto[1][0] ), &( zHisto[1][0] ),
665 delete allSlices[ sliceIndex ];
687 for (
long zpm=0; zpm<1 || (
m_chargeAware && zpm<2 ) ; ++zpm ) {
691 if((*nfiter).empty())
continue;
692 if((*nfiter).size() <= 2 )
continue;
698 pedestal += *(niter) + *(niter-1) + *(niter-2);
714 std::vector<double> zoutput;
715 std::vector<double> woutput;
723 long bending=0, bestPhi=0;
726 for (
long zpm=0; zpm<1 || (
m_chargeAware && zpm<2 ) ; ++zpm ) {
730 if((*nfiter).empty())
continue;
731 if((*nfiter).size() <= 2 )
continue;
735 ztest = *(niter-2) + *(niter-1) + *(niter);
736 if ( ztest <= 0 || ztest < maxh )
continue;
739 if ( ztest >= maxh ) {
740 long bintest = niter-(*nfiter).begin()-1;
743 (
m_preferCentralZ && std::abs( bintest - HalfZhistoBins ) < std::abs( binMax - HalfZhistoBins ) ) ) {
746 bestPhi = nfiter-
nHisto[zpm].begin();
762 double weightedMax = ( zHisto[bending][bestPhi][binMax] +
763 zHisto[bending][bestPhi][binMax+1] +
764 zHisto[bending][bestPhi][binMax-1] ) /maxh;
768 nHisto[bending][bestPhi][binMax] = -1;
769 nHisto[bending][bestPhi][binMax-1] = -1;
770 nHisto[bending][bestPhi][binMax+1] = -1;
771 zHisto[bending][bestPhi][binMax] = 0;
772 zHisto[bending][bestPhi][binMax-1] = 0;
773 zHisto[bending][bestPhi][binMax+1] = 0;
777 float dist2closestVtx = 1000;
778 for (
size_t iv = 0; iv < zoutput.size(); ++iv )
779 if ( fabs(weightedMax-zoutput[iv]) < dist2closestVtx ) {
781 dist2closestVtx = fabs(weightedMax-zoutput[iv]);
789 zoutput.push_back( weightedMax );
790 woutput.push_back( maxh );
824 for (
unsigned i=0 ;
i<zoutput.size() ;
i++ ) {
827 output->push_back( woutput[
i] - pedestal );
832 for (
unsigned i=0 ;
i<zoutput.size() ;
i++ )
output->push_back( zoutput[
i] );
842 for (
unsigned i=0 ;
i<zoutput.size() ;
i++ ) {
847 for (
unsigned i=0 ;
i<zoutput.size() ;
i++ ) {