DAClusterizerInZ_vect Class Reference

#include <DAClusterizerInZ_vect.h>

Inheritance diagram for DAClusterizerInZ_vect:

Classes
struct	track_t
struct	vertex_t

Public Member Functions
double	beta0 (const double betamax, track_t const &tks, vertex_t const &y) const
std::vector< std::vector < reco::TransientTrack > >	clusterize (const std::vector< reco::TransientTrack > &tracks) const
	DAClusterizerInZ_vect (const edm::ParameterSet &conf)
void	dump (const double beta, const vertex_t &y, const track_t &tks, const int verbosity=0) const
track_t	fill (const std::vector< reco::TransientTrack > &tracks) const
bool	merge (vertex_t &) const
bool	merge (vertex_t &y, double &beta) const
bool	purge (vertex_t &, track_t &, double &, const double) const
bool	split (const double beta, track_t &t, vertex_t &y) const
void	splitAll (vertex_t &y) const
double	update (double beta, track_t &gtracks, vertex_t &gvertices, bool useRho0, double &rho0) const
std::vector< TransientVertex >	vertices (const std::vector< reco::TransientTrack > &tracks, const int verbosity=0) const
Public Member Functions inherited from TrackClusterizerInZ
	TrackClusterizerInZ ()
	TrackClusterizerInZ (const edm::ParameterSet &conf)
virtual	~TrackClusterizerInZ ()

Private Attributes
float	betamax_
float	betastop_
double	coolingFactor_
double	d0CutOff_
double	dzCutOff_
int	maxIterations_
bool	useTc_
bool	verbose_
float	vertexSize_

Detailed Description

Description: separates event tracks into clusters along the beam line

Version which auto-vectorizes with gcc 4.6 or newer

Definition at line 21 of file DAClusterizerInZ_vect.h.

Constructor & Destructor Documentation

DAClusterizerInZ_vect::DAClusterizerInZ_vect

(

const edm::ParameterSet &

conf

)

Definition at line 15 of file DAClusterizerInZ_vect.cc.

References edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), LogDebug, and builder_last_value_cfg::Tmin.

                                                                        {

  // some defaults to avoid uninitialized variables
  verbose_ = conf.getUntrackedParameter<bool> ("verbose", false);
  
  betamax_ = 0.1;
  betastop_ = 1.0;
  /*
  coolingFactor_ = 0.6;
  maxIterations_ = 100;
  vertexSize_ = 0.01; // 0.1 mm
  dzCutOff_ = 4.0;  
  */
    // configure

  double Tmin = conf.getParameter<double> ("Tmin");
  vertexSize_ = conf.getParameter<double> ("vertexSize");
  coolingFactor_ = conf.getParameter<double> ("coolingFactor");
  useTc_=true;
  if(coolingFactor_<0){
    coolingFactor_=-coolingFactor_; 
    useTc_=false;
  }
  d0CutOff_ = conf.getParameter<double> ("d0CutOff");
  dzCutOff_ = conf.getParameter<double> ("dzCutOff");
  maxIterations_ = 100;
  if (Tmin == 0) {
    LogDebug("DAClusterizerinZ_vectorized")  << "DAClusterizerInZ: invalid Tmin" << Tmin
                         << "  reset do default " << 1. / betamax_ << endl;
  } else {
    betamax_ = 1. / Tmin;
  }

}

Member Function Documentation

double DAClusterizerInZ_vect::beta0	(	const double	betamax,
		track_t const &	tks,
		vertex_t const &	y
	)		const

Definition at line 355 of file DAClusterizerInZ_vect.cc.

References DAClusterizerInZ_vect::track_t::_dz2, DAClusterizerInZ_vect::track_t::_pi, DAClusterizerInZ_vect::track_t::_z, DAClusterizerInZ_vect::vertex_t::_z, a, b, DAClusterizerInZ_vect::track_t::GetSize(), DAClusterizerInZ_vect::vertex_t::GetSize(), i, relval_2017::k, dqm-mbProfile::log, nt, funct::pow(), and w.

                                                                                           {
  
  double T0 = 0; // max Tc for beta=0
  // estimate critical temperature from beta=0 (T=inf)
  const unsigned int nt = tks.GetSize();
  const unsigned int nv = y.GetSize();
  
  for (unsigned int k = 0; k < nv; k++) {
    
    // vertex fit at T=inf
    double sumwz = 0;
    double sumw = 0;
    for (unsigned int i = 0; i < nt; i++) {
      double w = tks._pi[i] * tks._dz2[i];
      sumwz += w * tks._z[i];
      sumw += w;
    }
    y._z[k] = sumwz / sumw;
    
    // estimate Tcrit, eventually do this in the same loop
    double a = 0, b = 0;
    for (unsigned int i = 0; i < nt; i++) {
      double dx = tks._z[i] - y._z[k];
      double w = tks._pi[i] * tks._dz2[i];
      a += w * std::pow(dx, 2) * tks._dz2[i];
      b += w;
    }
    double Tc = 2. * a / b; // the critical temperature of this vertex
    if (Tc > T0) T0 = Tc;
  }// vertex loop (normally there should be only one vertex at beta=0)
  
  if (T0 > 1. / betamax) {
    return betamax / std::pow(coolingFactor_, int(std::log(T0 * betamax) / std::log(coolingFactor_)) - 1);
  } else {
    // ensure at least one annealing step
    return betamax / coolingFactor_;
  }
}

vector< vector< reco::TransientTrack > > DAClusterizerInZ_vect::clusterize ( const std::vector< reco::TransientTrack > &  tracks ) const

virtual

Implements TrackClusterizerInZ.

Definition at line 660 of file DAClusterizerInZ_vect.cc.

References funct::abs(), HLT_25ns10e33_v2_cff::clusters, gather_cfg::cout, i, relval_2017::k, LogDebug, eostools::move(), position, MetAnalyzer::pv(), and HLT_25ns10e33_v2_cff::vertices.

                                                         {
  
  if (verbose_) {
    std::cout  << "###################################################" << endl;
    std::cout  << "# vectorized DAClusterizerInZ_vect::clusterize   nt=" << tracks.size() << endl;
    std::cout  << "###################################################" << endl;
  }
  
  vector<vector<reco::TransientTrack> > clusters;
  vector<TransientVertex> && pv = vertices(tracks);
  
  if (verbose_) {
    LogDebug("DAClusterizerinZ_vectorized")  << "# DAClusterizerInZ::clusterize   pv.size=" << pv.size()
                         << endl;
  }
  if (pv.size() == 0) {
    return clusters;
  }
  
  // fill into clusters and merge
  vector<reco::TransientTrack> aCluster = pv.begin()->originalTracks();
  
  for (auto k = pv.begin() + 1; k != pv.end(); k++) {
    if ( std::abs(k->position().z() - (k - 1)->position().z()) > (2 * vertexSize_)) {
      // close a cluster
      clusters.push_back(aCluster);
      aCluster.clear();
    }
    for (unsigned int i = 0; i < k->originalTracks().size(); i++) {
      aCluster.push_back(k->originalTracks()[i]);
    }
    
  }
  clusters.emplace_back(std::move(aCluster));
  
  return clusters;

}

void DAClusterizerInZ_vect::dump	(	const double	beta,
		const vertex_t &	y,
		const track_t &	tks,
		const int	verbosity = 0
	)		const

Definition at line 704 of file DAClusterizerInZ_vect.cc.

References DAClusterizerInZ_vect::track_t::_dz2, DAClusterizerInZ_vect::track_t::_pi, DAClusterizerInZ_vect::vertex_t::_pk, DAClusterizerInZ_vect::track_t::_z, DAClusterizerInZ_vect::vertex_t::_z, DAClusterizerInZ_vect::track_t::_Z_sum, beta, TauDecayModes::dec, Measurement1D::error(), create_public_lumi_plots::exp, F(), DAClusterizerInZ_vect::track_t::GetSize(), DAClusterizerInZ_vect::vertex_t::GetSize(), reco::TrackBase::highPurity, i, listHistos::IP, dqm-mbProfile::log, LogDebug, reco::HitPattern::MISSING_OUTER_HITS, nt, AlCaHLTBitMon_ParallelJobs::p, mathSSE::sqrt(), DAClusterizerInZ_vect::track_t::tt, and Measurement1D::value().

                                                  {

    const unsigned int nv = y.GetSize();
    const unsigned int nt = tks.GetSize();

    LogDebug("DAClusterizerinZ_vectorized")  << "-----DAClusterizerInZ::dump ----" << endl;
    LogDebug("DAClusterizerinZ_vectorized")  << "beta=" << beta << "   betamax= " << betamax_ << endl;
    LogDebug("DAClusterizerinZ_vectorized")  << "                                                               z= ";
    LogDebug("DAClusterizerinZ_vectorized")  << setprecision(4);
    for (unsigned int ivertex = 0; ivertex < nv; ++ ivertex) {
        LogDebug("DAClusterizerinZ_vectorized")  << setw(8) << fixed << y._z[ivertex];
    }
    LogDebug("DAClusterizerinZ_vectorized")  << endl << "T=" << setw(15) << 1. / beta
            << "                                            Tc= ";
    LogDebug("DAClusterizerinZ_vectorized")  << endl
            << "                                                               pk=";
    double sumpk = 0;
    for (unsigned int ivertex = 0; ivertex < nv; ++ ivertex) {
        LogDebug("DAClusterizerinZ_vectorized")  << setw(8) << setprecision(3) << fixed << y._pk[ivertex];
        sumpk += y._pk[ivertex];
    }
    LogDebug("DAClusterizerinZ_vectorized")  << endl;

    if (verbosity > 0) {
        double E = 0, F = 0;
        LogDebug("DAClusterizerinZ_vectorized")  << endl;
        LogDebug("DAClusterizerinZ_vectorized") 
        << "----       z +/- dz                ip +/-dip       pt    phi  eta    weights  ----"
        << endl;
        LogDebug("DAClusterizerinZ_vectorized")  << setprecision(4);
        for (unsigned int i = 0; i < nt; i++) {
            if (tks._Z_sum[i] > 0) {
                F -= std::log(tks._Z_sum[i]) / beta;
            }
            double tz = tks._z[i];
            LogDebug("DAClusterizerinZ_vectorized")  << setw(3) << i << ")" << setw(8) << fixed << setprecision(4)
             << tz << " +/-" << setw(6) << sqrt(1./tks._dz2[i]);

            if (tks.tt[i]->track().quality(reco::TrackBase::highPurity)) {
                LogDebug("DAClusterizerinZ_vectorized")  << " *";
            } else {
                LogDebug("DAClusterizerinZ_vectorized")  << "  ";
            }
            if (tks.tt[i]->track().hitPattern().hasValidHitInFirstPixelBarrel()) {
                LogDebug("DAClusterizerinZ_vectorized")  << "+";
            } else {
                LogDebug("DAClusterizerinZ_vectorized")  << "-";
            }
            LogDebug("DAClusterizerinZ_vectorized")  << setw(1)
             << tks.tt[i]->track().hitPattern().pixelBarrelLayersWithMeasurement(); // see DataFormats/TrackReco/interface/HitPattern.h
            LogDebug("DAClusterizerinZ_vectorized")  << setw(1)
             << tks.tt[i]->track().hitPattern().pixelEndcapLayersWithMeasurement();
            LogDebug("DAClusterizerinZ_vectorized")  << setw(1) << hex
                    << tks.tt[i]->track().hitPattern().trackerLayersWithMeasurement()
                    - tks.tt[i]->track().hitPattern().pixelLayersWithMeasurement()
                    << dec;
            LogDebug("DAClusterizerinZ_vectorized")  << "=" << setw(1) << hex
                    << tks.tt[i]->track().hitPattern().numberOfHits(reco::HitPattern::MISSING_OUTER_HITS)
                    << dec;

            Measurement1D IP =
                    tks.tt[i]->stateAtBeamLine().transverseImpactParameter();
            LogDebug("DAClusterizerinZ_vectorized")  << setw(8) << IP.value() << "+/-" << setw(6) << IP.error();
            LogDebug("DAClusterizerinZ_vectorized")  << " " << setw(6) << setprecision(2)
             << tks.tt[i]->track().pt() * tks.tt[i]->track().charge();
            LogDebug("DAClusterizerinZ_vectorized")  << " " << setw(5) << setprecision(2)
             << tks.tt[i]->track().phi() << " " << setw(5)
             << setprecision(2) << tks.tt[i]->track().eta();

            double sump = 0.;
            for (unsigned int ivertex = 0; ivertex < nv; ++ ivertex) {
                if ((tks._pi[i] > 0) && (tks._Z_sum[i] > 0)) {
                    //double p=pik(beta,tks[i],*k);
                    double p = y._pk[ivertex]  * exp(-beta * Eik(tks._z[i], y._z[ivertex], tks._dz2[i])) / tks._Z_sum[i];
                    if (p > 0.0001) {
                        LogDebug("DAClusterizerinZ_vectorized")  << setw(8) << setprecision(3) << p;
                    } else {
                        LogDebug("DAClusterizerinZ_vectorized")  << "    .   ";
                    }
                    E += p * Eik(tks._z[i], y._z[ivertex], tks._dz2[i]);
                    sump += p;
                } else {
                    LogDebug("DAClusterizerinZ_vectorized")  << "        ";
                }
            }
            LogDebug("DAClusterizerinZ_vectorized")  << endl;
        }
        LogDebug("DAClusterizerinZ_vectorized")  << endl << "T=" << 1 / beta << " E=" << E << " n=" << y.GetSize()
             << "  F= " << F << endl << "----------" << endl;
    }
}

DAClusterizerInZ_vect::track_t DAClusterizerInZ_vect::fill

(

const std::vector< reco::TransientTrack > &

tracks

)

const

Definition at line 63 of file DAClusterizerInZ_vect.cc.

References funct::abs(), DAClusterizerInZ_vect::track_t::AddItem(), SiPixelRawToDigiRegional_cfi::beamSpot, HLT_25ns10e33_v2_cff::beamspot, reco::BeamSpot::BeamWidthX(), reco::BeamSpot::BeamWidthY(), geometryDiff::epsilon, Measurement1D::error(), DAClusterizerInZ_vect::track_t::ExtractRaw(), DAClusterizerInZ_vect::track_t::GetSize(), edm::isNotFinite(), LogDebug, LogTrace, min(), position, funct::pow(), and Measurement1D::value().

                                                                           {

  // prepare track data for clustering
  LogDebug("DAClusterizerinZ_vectorized") << "input values"; 
  track_t tks;
  for (auto it = tracks.begin(); it!= tracks.end(); it++){
    if (!(*it).isValid()) continue;
    double t_pi=1.;
    double t_z = ((*it).stateAtBeamLine().trackStateAtPCA()).position().z();
    if (std::abs(t_z) > 1000.) continue;
    auto const & t_mom = (*it).stateAtBeamLine().trackStateAtPCA().momentum();
    //  get the beam-spot
    reco::BeamSpot beamspot = (it->stateAtBeamLine()).beamSpot();
    double t_dz2 = 
      std::pow((*it).track().dzError(), 2) // track errror
      + (std::pow(beamspot.BeamWidthX()*t_mom.x(),2)+std::pow(beamspot.BeamWidthY()*t_mom.y(),2))*std::pow(t_mom.z(),2)/std::pow(t_mom.perp2(),2) // beam spot width
      + std::pow(vertexSize_, 2); // intrinsic vertex size, safer for outliers and short lived decays
    t_dz2 = 1./ t_dz2;
    if (edm::isNotFinite(t_dz2) || t_dz2 < std::numeric_limits<double>::min() ) continue;
    if (d0CutOff_ > 0) {
      Measurement1D atIP =
    (*it).stateAtBeamLine().transverseImpactParameter();// error constains beamspot
      t_pi = 1. / (1. + local_exp(std::pow(atIP.value() / atIP.error(), 2) - std::pow(d0CutOff_, 2))); // reduce weight for high ip tracks
      if (edm::isNotFinite(t_pi) ||  t_pi < std::numeric_limits<double>::epsilon())  continue; // usually is > 0.99
    }
    LogTrace("DAClusterizerinZ_vectorized") << t_z <<' '<< t_dz2 <<' '<< t_pi;
    tks.AddItem(t_z, t_dz2, &(*it), t_pi);
  }
  tks.ExtractRaw();
  
  if (verbose_) {
    LogDebug("DAClusterizerinZ_vectorized") << "Track count " << tks.GetSize() << std::endl;
  }
  
  return tks;
}

bool DAClusterizerInZ_vect::merge

(

vertex_t &

y

)

const

Definition at line 278 of file DAClusterizerInZ_vect.cc.

References DAClusterizerInZ_vect::vertex_t::_pk, DAClusterizerInZ_vect::vertex_t::_z, funct::abs(), DAClusterizerInZ_vect::vertex_t::GetSize(), relval_2017::k, LogDebug, and DAClusterizerInZ_vect::vertex_t::RemoveItem().

                                                    {
  // merge clusters that collapsed or never separated, return true if vertices were merged, false otherwise
  
  const unsigned int nv = y.GetSize();
  
  if (nv < 2) return false;

  for (unsigned int k = 0; (k + 1) < nv; k++) {
    //if ((k+1)->z - k->z<1.e-2){  // note, no fabs here, maintains z-ordering  (with split()+merge() at every temperature)
    if (std::abs(y._z[k + 1] - y._z[k]) < 1.e-2) { // with fabs if only called after freeze-out (splitAll() at highter T)
      y._pk[k] += y._pk[k + 1];
      y._z[k] = 0.5 * (y._z[k] + y._z[k + 1]);
      
      y.RemoveItem(k + 1);
      
      if ( verbose_)
    {
      LogDebug("DAClusterizerinZ_vectorized") << "Merging vertices k = " << k << std::endl;
    }
      
      return true;
    }
  }
  
  return false;
}

bool DAClusterizerInZ_vect::merge	(	vertex_t &	y,
		double &	beta
	)		const

Definition at line 240 of file DAClusterizerInZ_vect.cc.

References DAClusterizerInZ_vect::vertex_t::_pk, DAClusterizerInZ_vect::vertex_t::_sw, DAClusterizerInZ_vect::vertex_t::_swE, DAClusterizerInZ_vect::vertex_t::_z, DAClusterizerInZ_vect::vertex_t::GetSize(), relval_2017::k, funct::pow(), DAClusterizerInZ_vect::vertex_t::RemoveItem(), and rho.

                                                                 {
  // merge clusters that collapsed or never separated,
  // only merge if the estimated critical temperature of the merged vertex is below the current temperature
  // return true if vertices were merged, false otherwise
  const unsigned int nv = y.GetSize();

  if (nv < 2)
    return false;


  for (unsigned int k = 0; (k + 1) < nv; k++) {
    if (fabs(y._z[k + 1] - y._z[k]) < 2.e-2) {
      double rho=y._pk[k] + y._pk[k+1];
      double swE=y._swE[k]+y._swE[k+1] - y._pk[k]*y._pk[k+1] /rho *std::pow(y._z[k+1]-y._z[k],2);
      double Tc=2*swE/(y._sw[k]+y._sw[k+1]);

      if(Tc*beta<1){
        if(rho>0){
      y._z[k] = (y._pk[k]*y._z[k] + y._pk[k+1]*y._z[k + 1])/rho;
        }else{
      y._z[k] = 0.5 * (y._z[k] + y._z[k + 1]);
        }
        y._pk[k] = rho;
        y._sw[k]+=y._sw[k+1];
        y._swE[k]=swE;
        y.RemoveItem(k+1);
        return true;
      }
    }
  }

  return false;
}

bool DAClusterizerInZ_vect::purge	(	vertex_t &	y,
		track_t &	tks,
		double &	rho0,
		const double	beta
	)		const

Definition at line 306 of file DAClusterizerInZ_vect.cc.

References DAClusterizerInZ_vect::track_t::_dz2, DAClusterizerInZ_vect::track_t::_pi, DAClusterizerInZ_vect::vertex_t::_pk, DAClusterizerInZ_vect::track_t::_z, DAClusterizerInZ_vect::vertex_t::_z, DAClusterizerInZ_vect::track_t::_Z_sum, DAClusterizerInZ_vect::track_t::GetSize(), DAClusterizerInZ_vect::vertex_t::GetSize(), i, relval_2017::k, reco::ParticleMasses::k0, LogDebug, nt, AlCaHLTBitMon_ParallelJobs::p, and DAClusterizerInZ_vect::vertex_t::RemoveItem().

                                                                                                {
  // eliminate clusters with only one significant/unique track
  const unsigned int nv = y.GetSize();
  const unsigned int nt = tks.GetSize();

  if (nv < 2)
    return false;
  
  double sumpmin = nt;
  unsigned int k0 = nv;
  
  for (unsigned int k = 0; k < nv; k++) {
    
    int nUnique = 0;
    double sump = 0;
    double pmax = y._pk[k] / (y._pk[k] + rho0 * local_exp(-beta * dzCutOff_* dzCutOff_));
    
    for (unsigned int i = 0; i < nt; i++) {
      
      if (tks._Z_sum[i] > 0) {
    //double p=pik(beta,tks[i],*k);
    double p = y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k], tks._dz2[i])) / tks._Z_sum[i];
    sump += p;
    if ((p > 0.9 * pmax) && (tks._pi[i] > 0)) {
      nUnique++;
    }
            }
    }
    
    if ((nUnique < 2) && (sump < sumpmin)) {
      sumpmin = sump;
      k0 = k;
    }
  }
  
  if (k0 != nv) {
    if (verbose_) {
      LogDebug("DAClusterizerinZ_vectorized")  << "eliminating prototype at " << y._z[k0] << " with sump="
                           << sumpmin << endl;
    }
    //rho0+=k0->pk;
    y.RemoveItem(k0);
    return true;
  } else {
    return false;
  }
}

bool DAClusterizerInZ_vect::split	(	const double	beta,
		track_t &	t,
		vertex_t &	y
	)		const

Definition at line 396 of file DAClusterizerInZ_vect.cc.

References DAClusterizerInZ_vect::track_t::_dz2, DAClusterizerInZ_vect::vertex_t::_pk, DAClusterizerInZ_vect::vertex_t::_sw, DAClusterizerInZ_vect::vertex_t::_swE, DAClusterizerInZ_vect::track_t::_z, DAClusterizerInZ_vect::vertex_t::_z, DAClusterizerInZ_vect::track_t::_Z_sum, hlt_jetmet_dqm_QT_fromfile_cfg::critical, alignCSCRings::e, geometryDiff::epsilon, DAClusterizerInZ_vect::track_t::GetSize(), DAClusterizerInZ_vect::vertex_t::GetSize(), i, DAClusterizerInZ_vect::vertex_t::InsertItem(), relval_2017::k, nt, AlCaHLTBitMon_ParallelJobs::p, p1, p2, edm::second(), split, w, and w2.

                                                                                 {
  // split only critical vertices (Tc >~ T=1/beta   <==>   beta*Tc>~1)
  // an update must have been made just before doing this (same beta, no merging)
  // returns true if at least one cluster was split
  
  double epsilon=1e-3;      // split all single vertices by 10 um
  unsigned int nv = y.GetSize();
  
  // avoid left-right biases by splitting highest Tc first
  
  std::vector<std::pair<double, unsigned int> > critical;
  for(unsigned int k=0; k<nv; k++){
    double Tc= 2*y._swE[k]/y._sw[k];
    if (beta*Tc > 1.){
      critical.push_back( make_pair(Tc, k));
    }
  }
  if (critical.size()==0) return false;
  std::stable_sort(critical.begin(), critical.end(), std::greater<std::pair<double, unsigned int> >() );
  
  
  bool split=false;
  const unsigned int nt = tks.GetSize();

  for(unsigned int ic=0; ic<critical.size(); ic++){
    unsigned int k=critical[ic].second;
    // estimate subcluster positions and weight
    double p1=0, z1=0, w1=0;
    double p2=0, z2=0, w2=0;
    for(unsigned int i=0; i<nt; i++){
      if (tks._Z_sum[i] > 0) {
    double p = y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k], tks._dz2[i])) / tks._Z_sum[i];
        double w=p*tks._dz2[i];
        if(tks._z[i]<y._z[k]){
          p1+=p; z1+=w*tks._z[i]; w1+=w;
        }else{
          p2+=p; z2+=w*tks._z[i]; w2+=w;
        }
      }
    }
    if(w1>0){  z1=z1/w1;} else{z1=y._z[k]-epsilon;}
    if(w2>0){  z2=z2/w2;} else{z2=y._z[k]+epsilon;}

    // reduce split size if there is not enough room
    if( ( k   > 0 ) && ( y._z[k-1]>=z1 ) ){ z1=0.5*(y._z[k]+y._z[k-1]); }
    if( ( k+1 < nv) && ( y._z[k+1]<=z2 ) ){ z2=0.5*(y._z[k]+y._z[k+1]); }

    // split if the new subclusters are significantly separated
    if( (z2-z1)>epsilon){
      split=true;
      double pk1=p1*y._pk[k]/(p1+p2);
      double pk2=p2*y._pk[k]/(p1+p2);
      y._z[k]  =  z2;
      y._pk[k] = pk2;
      y.InsertItem(k, z1, pk1);
      nv++;

     // adjust remaining pointers
      for(unsigned int jc=ic; jc<critical.size(); jc++){
        if (critical[jc].second>k) {critical[jc].second++;}
      }
    }
  }
  return split;
}

void DAClusterizerInZ_vect::splitAll

(

vertex_t &

y

)

const

Definition at line 464 of file DAClusterizerInZ_vect.cc.

References DAClusterizerInZ_vect::vertex_t::_pk, DAClusterizerInZ_vect::vertex_t::_z, DAClusterizerInZ_vect::vertex_t::AddItem(), DAClusterizerInZ_vect::vertex_t::DebugOut(), alignCSCRings::e, geometryDiff::epsilon, DAClusterizerInZ_vect::vertex_t::ExtractRaw(), DAClusterizerInZ_vect::vertex_t::GetSize(), relval_2017::k, LogDebug, DAClusterizerInZ_vect::vertex_t::pk, and DAClusterizerInZ_vect::vertex_t::z.

                                                        {

  const unsigned int nv = y.GetSize();
  
  double epsilon = 1e-3; // split all single vertices by 10 um
  double zsep = 2 * epsilon; // split vertices that are isolated by at least zsep (vertices that haven't collapsed)
  vertex_t y1;
  
  if (verbose_) {
    LogDebug("DAClusterizerinZ_vectorized") << "Before Split "<< std::endl;
    y.DebugOut();
  }

  for (unsigned int k = 0; k < nv; k++) {
    if (
    ( (k == 0)          || ( y._z[k - 1]    < (y._z[k] - zsep)) ) &&
    ( ((k + 1) == nv)   || ( y._z[k + 1]    > (y._z[k] + zsep)) )   )
      {
    // isolated prototype, split
    double new_z = y.z[k] - epsilon;
    y.z[k] = y.z[k] + epsilon;
    
    double new_pk = 0.5 * y.pk[k];
    y.pk[k] = 0.5 * y.pk[k];
    
    y1.AddItem(new_z, new_pk);
    y1.AddItem(y._z[k], y._pk[k]);  
      }
    else if ( (y1.GetSize() == 0 ) ||
          (y1._z[y1.GetSize() - 1] <  (y._z[k] - zsep)  ))
      {
    y1.AddItem(y._z[k], y._pk[k]);
      }
    else
      {
    y1._z[y1.GetSize() - 1] = y1._z[y1.GetSize() - 1] - epsilon;
    y._z[k] = y._z[k] + epsilon;
    y1.AddItem( y._z[k] , y._pk[k]);
      }
  }// vertex loop

  y = y1;
  y.ExtractRaw();
  
  if (verbose_) {
    LogDebug("DAClusterizerinZ_vectorized") << "After split " << std::endl;
    y.DebugOut();
  }
}

double DAClusterizerInZ_vect::update	(	double	beta,
		track_t &	gtracks,
		vertex_t &	gvertices,
		bool	useRho0,
		double &	rho0
	)		const

Definition at line 108 of file DAClusterizerInZ_vect.cc.

References DAClusterizerInZ_vect::vertex_t::_ei, DAClusterizerInZ_vect::vertex_t::_ei_cache, DAClusterizerInZ_vect::track_t::_pi, DAClusterizerInZ_vect::vertex_t::_se, DAClusterizerInZ_vect::vertex_t::_sw, DAClusterizerInZ_vect::vertex_t::_swE, DAClusterizerInZ_vect::vertex_t::_swz, DAClusterizerInZ_vect::track_t::_Z_sum, beta, delta, DAClusterizerInZ_vect::track_t::GetSize(), DAClusterizerInZ_vect::vertex_t::GetSize(), relval_2017::k, LogDebug, nt, funct::pow(), testEve_cfg::tracks, HLT_25ns10e33_v2_cff::vertices, and w.

Referenced by progressbar.ProgressBar::__next__(), MatrixUtil.Matrix::__setitem__(), MatrixUtil.Steps::__setitem__(), Vispa.Gui.VispaWidget.VispaWidget::autosize(), Vispa.Views.LineDecayView.LineDecayContainer::createObject(), Vispa.Views.LineDecayView.LineDecayContainer::deselectAllObjects(), Vispa.Gui.VispaWidgetOwner.VispaWidgetOwner::deselectAllWidgets(), Vispa.Gui.VispaWidget.VispaWidget::enableAutosizing(), progressbar.ProgressBar::finish(), Vispa.Gui.MenuWidget.MenuWidget::leaveEvent(), Vispa.Gui.VispaWidgetOwner.VispaWidgetOwner::mouseMoveEvent(), Vispa.Gui.MenuWidget.MenuWidget::mouseMoveEvent(), Vispa.Views.LineDecayView.LineDecayContainer::mouseMoveEvent(), Vispa.Gui.VispaWidgetOwner.VispaWidgetOwner::mouseReleaseEvent(), Vispa.Views.LineDecayView.LineDecayContainer::objectMoved(), MatrixUtil.Steps::overwrite(), Vispa.Views.LineDecayView.LineDecayContainer::removeObject(), Vispa.Gui.ConnectableWidget.ConnectableWidget::removePorts(), Vispa.Gui.FindDialog.FindDialog::reset(), Vispa.Gui.PortConnection.PointToPointConnection::select(), Vispa.Gui.VispaWidget.VispaWidget::select(), Vispa.Views.LineDecayView.LineDecayContainer::select(), Vispa.Gui.VispaWidget.VispaWidget::setText(), Vispa.Gui.VispaWidget.VispaWidget::setTitle(), Vispa.Gui.ZoomableWidget.ZoomableWidget::setZoom(), Vispa.Views.LineDecayView.LineDecayContainer::setZoom(), and Vispa.Gui.PortConnection.PointToPointConnection::updateConnection().

                                                                              {

  //update weights and vertex positions
  // mass constrained annealing without noise
  // returns the squared sum of changes of vertex positions
  
  const unsigned int nt = gtracks.GetSize();
  const unsigned int nv = gvertices.GetSize();
  
  //initialize sums
  double sumpi = 0;
  
  // to return how much the prototype moved
  double delta = 0;
  

  // intial value of a sum
  double Z_init = 0;
  
  // define kernels
  auto kernel_calc_exp_arg = [ beta, nv ] ( const unsigned int itrack,
                         track_t const& tracks,
                         vertex_t const& vertices ) {
    const double track_z = tracks._z[itrack];
    const double botrack_dz2 = -beta*tracks._dz2[itrack];

    // auto-vectorized
    for ( unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
      auto mult_res =  track_z - vertices._z[ivertex];
      vertices._ei_cache[ivertex] = botrack_dz2 * ( mult_res * mult_res );
    }
  };
  
  auto kernel_add_Z = [ nv, Z_init ] (vertex_t const& vertices) -> double
    {
      double ZTemp = Z_init;
      for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) { 
    ZTemp += vertices._pk[ivertex] * vertices._ei[ivertex];
      }
      return ZTemp;
    };

  auto kernel_calc_normalization = [ beta, nv ] (const unsigned int track_num,
                          track_t & tks_vec,
                          vertex_t & y_vec ) {
    auto tmp_trk_pi = tks_vec._pi[track_num];
    auto o_trk_Z_sum = 1./tks_vec._Z_sum[track_num];
    auto o_trk_dz2 = tks_vec._dz2[track_num];
    auto tmp_trk_z = tks_vec._z[track_num];
    auto obeta =  -1./beta;
    
    // auto-vectorized
    for (unsigned int k = 0; k < nv; ++k) {
      y_vec._se[k] +=  y_vec._ei[k] * (tmp_trk_pi* o_trk_Z_sum);
      auto w = y_vec._pk[k] * y_vec._ei[k] * (tmp_trk_pi*o_trk_Z_sum *o_trk_dz2);
      y_vec._sw[k]  += w;
      y_vec._swz[k] += w * tmp_trk_z;
      y_vec._swE[k] += w * y_vec._ei_cache[k]*obeta;
    }
  };
  
  
  for (auto ivertex = 0U; ivertex < nv; ++ivertex) {
    gvertices._se[ivertex] = 0.0;
    gvertices._sw[ivertex] = 0.0;
    gvertices._swz[ivertex] = 0.0;
    gvertices._swE[ivertex] = 0.0;
  }
  
  
  // independpent of loop
  if ( useRho0 )
    {
      Z_init = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_); // cut-off
    }
  
  // loop over tracks
  for (auto itrack = 0U; itrack < nt; ++itrack) {
    kernel_calc_exp_arg(itrack, gtracks, gvertices);
    local_exp_list(gvertices._ei_cache, gvertices._ei, nv);
    
    gtracks._Z_sum[itrack] = kernel_add_Z(gvertices);
    
    // used in the next major loop to follow
    if (!useRho0)
      sumpi += gtracks._pi[itrack];
    
    if (gtracks._Z_sum[itrack] > 0) {
      kernel_calc_normalization(itrack, gtracks, gvertices);
    }
  }
  
  // now update z and pk
  auto kernel_calc_z = [  sumpi, nv, this, useRho0 ] (vertex_t & vertices ) -> double {
    
    double delta=0;
    // does not vectorizes
    for (unsigned int ivertex = 0; ivertex < nv; ++ ivertex ) {
      if (vertices._sw[ivertex] > 0) {
    auto znew = vertices._swz[ ivertex ] / vertices._sw[ ivertex ];
    // prevents from vectorizing if 
    delta += std::pow( vertices._z[ ivertex ] - znew, 2 );
    vertices._z[ ivertex ] = znew;
      }
#ifdef VI_DEBUG
      else {
    edm::LogInfo("sumw") << "invalid sum of weights in fit: " << vertices._sw[ivertex] << endl;
    if (this->verbose_) {
      LogDebug("DAClusterizerinZ_vectorized")  << " a cluster melted away ?  pk=" << vertices._pk[ ivertex ] << " sumw="
                           << vertices._sw[ivertex] << endl;
    }
      }
#endif
    }
    // dont do, if rho cut
    if ( ! useRho0 ){
      auto osumpi = 1./sumpi;
      for (unsigned int ivertex = 0; ivertex < nv; ++ ivertex )
    vertices._pk[ ivertex ] = vertices._pk[ ivertex ] * vertices._se[ ivertex ] * osumpi;
    }
    return delta;
  };
  
  delta += kernel_calc_z(gvertices);
  
  // return how much the prototypes moved
  return delta;
}

vector< TransientVertex > DAClusterizerInZ_vect::vertices	(	const std::vector< reco::TransientTrack > &	tracks,
		const int	verbosity = 0
	)		const

Definition at line 517 of file DAClusterizerInZ_vect.cc.

References DAClusterizerInZ_vect::track_t::_dz2, DAClusterizerInZ_vect::track_t::_pi, DAClusterizerInZ_vect::vertex_t::_pk, DAClusterizerInZ_vect::track_t::_z, DAClusterizerInZ_vect::vertex_t::_z, DAClusterizerInZ_vect::track_t::_Z_sum, DAClusterizerInZ_vect::vertex_t::AddItem(), beta, HLT_25ns10e33_v2_cff::clusters, hcal_timing_source_file_cfg::dump, alignCSCRings::e, DAClusterizerInZ_vect::track_t::ExtractRaw(), lumiContext::fill, DAClusterizerInZ_vect::track_t::GetSize(), DAClusterizerInZ_vect::vertex_t::GetSize(), i, edm::isNotFinite(), relval_2017::k, LogDebug, MatrixUtil::merge(), nt, AlCaHLTBitMon_ParallelJobs::p, split, DAClusterizerInZ_vect::track_t::tt, update, and findQualityFiles::v.

                                                                                                    {
  track_t && tks = fill(tracks);
  tks.ExtractRaw();
  
  unsigned int nt = tks.GetSize();
  double rho0 = 0.0; // start with no outlier rejection
  
  vector<TransientVertex> clusters;
  if (tks.GetSize() == 0) return clusters;
  
  vertex_t y; // the vertex prototypes
  
  // initialize:single vertex at infinite temperature
  y.AddItem( 0, 1.0);
  
  int niter = 0; // number of iterations
  
  
  // estimate first critical temperature
  double beta = beta0(betamax_, tks, y);
  if ( verbose_) LogDebug("DAClusterizerinZ_vectorized") << "Beta0 is " << beta << std::endl;
  
  niter = 0;
  while ((update(beta, tks, y, false, rho0) > 1.e-6) && 
     (niter++ < maxIterations_)) {}

  // annealing loop, stop when T<Tmin  (i.e. beta>1/Tmin)
  while (beta < betamax_) {
    if(useTc_){
      update(beta, tks,y, false, rho0);
      while(merge(y,beta)){update(beta, tks,y, false, rho0);}
      split(beta, tks,y);
      beta=beta/coolingFactor_;
    }else{
      beta=beta/coolingFactor_;
      splitAll(y);
    }
    
    // make sure we are not too far from equilibrium before cooling further
    niter = 0;
    while ((update(beta, tks, y, false, rho0) > 1.e-6) && 
       (niter++ < maxIterations_)) {}
  }
  

  if(useTc_){
    // last round of splitting, make sure no critical clusters are left
    update(beta, tks,y, false, rho0);// make sure Tc is up-to-date
    while(merge(y,beta)){update(beta, tks,y, false, rho0);}
    unsigned int ntry=0;
    while( split(beta, tks,y) && (ntry++<10) ){
      niter=0; 
      while((update(beta, tks,y, false, rho0)>1.e-6)  && (niter++ < maxIterations_)){}
      merge(y,beta);
      update(beta, tks,y, false, rho0);
    }
  }else{
    // merge collapsed clusters 
    while(merge(y,beta)){update(beta, tks,y, false, rho0);}  
    //while(merge(y)){}   original code
  }
  
  if (verbose_) {
    LogDebug("DAClusterizerinZ_vectorized")  << "dump after 1st merging " << endl;
    dump(beta, y, tks, 2);
  }
  
  // switch on outlier rejection
  rho0 = 1. / nt;
  
  // auto-vectorized
  for (unsigned int k = 0; k < y.GetSize(); k++) y._pk[k] = 1.; // democratic
  niter = 0;
  while ((update(beta, tks, y, true, rho0) > 1.e-8) && (niter++ < maxIterations_)) {}
  if (verbose_) {
    LogDebug("DAClusterizerinZ_vectorized")  << "rho0=" << rho0 << " niter=" << niter << endl;
    dump(beta, y, tks, 2);
  }
  
  // merge again  (some cluster split by outliers collapse here)
  while (merge(y)) {}
  if (verbose_) {
    LogDebug("DAClusterizerinZ_vectorized")  << "dump after 2nd merging " << endl;
    dump(beta, y, tks, 2);
  }

  // continue from freeze-out to Tstop (=1) without splitting, eliminate insignificant vertices
  while (beta <= betastop_) {
    while (purge(y, tks, rho0, beta)) {
      niter = 0;
      while ((update(beta, tks, y, true, rho0) > 1.e-6) && (niter++ < maxIterations_)) {}
    }
    beta /= coolingFactor_;
    niter = 0;
    while ((update(beta, tks, y, true, rho0) > 1.e-6) && (niter++< maxIterations_)) {}
  }

  if (verbose_) {
    LogDebug("DAClusterizerinZ_vectorized")  << "Final result, rho0=" << rho0 << endl;
    dump(beta, y, tks, 2);
  }

  // select significant tracks and use a TransientVertex as a container
  GlobalError dummyError;
  
  // ensure correct normalization of probabilities, should makes double assginment reasonably impossible
  const unsigned int nv = y.GetSize();
  for (unsigned int k = 0; k < nv; k++)
     if ( edm::isNotFinite(y._pk[k]) || edm::isNotFinite(y._z[k]) ) { y._pk[k]=0; y._z[k]=0;}

  for (unsigned int i = 0; i < nt; i++) // initialize
    tks._Z_sum[i] = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_);

  // improve vectorization (does not require reduction ....)
  for (unsigned int k = 0; k < nv; k++) {
     for (unsigned int i = 0; i < nt; i++)  
      tks._Z_sum[i] += y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k],tks._dz2[i]));
  }


  for (unsigned int k = 0; k < nv; k++) {
    GlobalPoint pos(0, 0, y._z[k]);
    
    vector<reco::TransientTrack> vertexTracks;
    for (unsigned int i = 0; i < nt; i++) {
      if (tks._Z_sum[i] > 0) {
    
    double p = y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k],
                            tks._dz2[i])) / tks._Z_sum[i];
    if ((tks._pi[i] > 0) && (p > 0.5)) {
      vertexTracks.push_back(*(tks.tt[i]));
      tks._Z_sum[i] = 0;
    } // setting Z=0 excludes double assignment
      }
    }
    TransientVertex v(pos, dummyError, vertexTracks, 0);
    clusters.push_back(v);
  }

  return clusters;

}

Member Data Documentation

float DAClusterizerInZ_vect::betamax_

private

Definition at line 202 of file DAClusterizerInZ_vect.h.

float DAClusterizerInZ_vect::betastop_

private

Definition at line 203 of file DAClusterizerInZ_vect.h.

double DAClusterizerInZ_vect::coolingFactor_

private

Definition at line 201 of file DAClusterizerInZ_vect.h.

double DAClusterizerInZ_vect::d0CutOff_

private

Definition at line 205 of file DAClusterizerInZ_vect.h.

double DAClusterizerInZ_vect::dzCutOff_

private

Definition at line 204 of file DAClusterizerInZ_vect.h.

int DAClusterizerInZ_vect::maxIterations_

private

Definition at line 200 of file DAClusterizerInZ_vect.h.

bool DAClusterizerInZ_vect::useTc_

private

Definition at line 206 of file DAClusterizerInZ_vect.h.

bool DAClusterizerInZ_vect::verbose_

private

Definition at line 198 of file DAClusterizerInZ_vect.h.

float DAClusterizerInZ_vect::vertexSize_

private

Definition at line 199 of file DAClusterizerInZ_vect.h.