DAClusterizerInZ Class Reference

#include <DAClusterizerInZ.h>

Inheritance diagram for DAClusterizerInZ:

Classes
struct	track_t
struct	vertex_t

Public Member Functions
double	beta0 (const double betamax, std::vector< track_t > &tks, std::vector< vertex_t > &y) const
std::vector< std::vector < reco::TransientTrack > >	clusterize (const std::vector< reco::TransientTrack > &tracks) const
	DAClusterizerInZ (const edm::ParameterSet &conf)
void	dump (const double beta, const std::vector< vertex_t > &y, const std::vector< track_t > &tks, const int verbosity=0) const
double	Eik (const track_t &t, const vertex_t &k) const
std::vector< track_t >	fill (const std::vector< reco::TransientTrack > &tracks) const
bool	merge (std::vector< vertex_t > &, int) const
bool	purge (std::vector< vertex_t > &, std::vector< track_t > &, double &, const double) const
void	splitAll (std::vector< track_t > &tks, std::vector< vertex_t > &y) const
double	update (double beta, std::vector< track_t > &tks, std::vector< vertex_t > &y) const
double	update (double beta, std::vector< track_t > &tks, std::vector< vertex_t > &y, double &) 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	verbose_
float	vertexSize_

Detailed Description

Description: separates event tracks into clusters along the beam line

Definition at line 21 of file DAClusterizerInZ.h.

Constructor & Destructor Documentation

DAClusterizerInZ::DAClusterizerInZ

(

const edm::ParameterSet &

conf

)

Definition at line 348 of file DAClusterizerInZ.cc.

References gather_cfg::cout, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), 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.8;
  maxIterations_=100;
  vertexSize_=0.05;  // 0.5 mm
  dzCutOff_=4.0;   // Adaptive Fitter uses 3.0 but that appears to be a bit tight here sometimes

  // configure

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

}

Member Function Documentation

double DAClusterizerInZ::beta0	(	const double	betamax,
		std::vector< track_t > &	tks,
		std::vector< vertex_t > &	y
	)		const

Definition at line 269 of file DAClusterizerInZ.cc.

References a, b, i, gen::k, funct::log(), nt, and funct::pow().

                         {
  
  double T0=0;  // max Tc for beta=0
  // estimate critical temperature from beta=0 (T=inf)
  unsigned int nt=tks.size();

  for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){

    // vertex fit at T=inf 
    double sumwz=0;
    double sumw=0;
    for(unsigned int i=0; i<nt; i++){
      double w=tks[i].pi/tks[i].dz2;
      sumwz+=w*tks[i].z;
      sumw +=w;
    }
    k->z=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[i].z-(k->z);
      double w=tks[i].pi/tks[i].dz2;
      a+=w*pow(dx,2)/tks[i].dz2;
      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/pow(coolingFactor_, int(log(T0*betamax)/log(coolingFactor_))-1 );
  }else{
    // ensure at least one annealing step
    return betamax/coolingFactor_;
  }
}

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

virtual

Implements TrackClusterizerInZ.

Definition at line 556 of file DAClusterizerInZ.cc.

References gather_cfg::cout, i, gen::k, and position.

{
  if(verbose_) {
    cout << "###################################################" << endl;
    cout << "# DAClusterizerInZ::clusterize   nt="<<tracks.size() << endl;
    cout << "###################################################" << endl;
  }

  vector< vector<reco::TransientTrack> > clusters;
  vector< TransientVertex > pv=vertices(tracks);

  if(verbose_){ cout << "# 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(vector<TransientVertex>::iterator k=pv.begin()+1; k!=pv.end(); k++){
    if ( 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().at(i)); }
    
  }
  clusters.push_back(aCluster);
  
  
  return clusters;

}

void DAClusterizerInZ::dump	(	const double	beta,
		const std::vector< vertex_t > &	y,
		const std::vector< track_t > &	tks,
		const int	verbosity = 0
	)		const

Definition at line 376 of file DAClusterizerInZ.cc.

References beta, gather_cfg::cout, Measurement1D::error(), funct::exp(), reco::TrackBase::highPurity, i, gen::k, funct::log(), L1TEmulatorMonitor_cff::p, pi, mathSSE::sqrt(), matplotRender::t, Measurement1D::value(), and Gflash::Z.

                                                                                                                          {

  // copy and sort for nicer printout
  vector<track_t> tks; 
  for(vector<track_t>::const_iterator t=tks0.begin(); t!=tks0.end(); t++){tks.push_back(*t); }
  stable_sort(tks.begin(), tks.end(), recTrackLessZ1);

  cout << "-----DAClusterizerInZ::dump ----" << endl;
  cout << "beta=" << beta << "   betamax= " << betamax_ << endl;
  cout << "                                                               z= ";
  cout.precision(4);
  for(vector<vertex_t>::const_iterator k=y.begin(); k!=y.end(); k++){
    cout  <<  setw(8) << fixed << k->z ;
  }
  cout << endl << "T=" << setw(15) << 1./beta <<"                                            Tc= ";
  cout << endl << "                                                               pk=";
  double sumpk=0;
  for(vector<vertex_t>::const_iterator k=y.begin(); k!=y.end(); k++){
    cout <<  setw(8) <<  setprecision(3) <<  fixed << k->pk;
    sumpk+=k->pk;
  }
  cout  << endl;

  if(verbosity>0){
    double E=0, F=0;
    cout << endl;
    cout << "----       z +/- dz                ip +/-dip       pt    phi  eta    weights  ----" << endl;
    cout.precision(4);
    for(unsigned int i=0; i<tks.size(); i++){
      if (tks[i].Z>0){  F-=log(tks[i].Z)/beta;}
      double tz= tks[i].z;
      cout <<  setw (3)<< i << ")" <<  setw (8) << fixed << setprecision(4)<<  tz << " +/-" <<  setw (6)<< sqrt(tks[i].dz2);

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

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

      double sump=0.;
      for(vector<vertex_t>::const_iterator k=y.begin(); k!=y.end(); k++){
    if((tks[i].pi>0)&&(tks[i].Z>0)){
      //double p=pik(beta,tks[i],*k);
      double p=k->pk * exp(-beta*Eik(tks[i],*k)) / tks[i].Z; 
      if( p > 0.0001){
        cout <<  setw (8) <<  setprecision(3) << p;
      }else{
        cout << "    .   ";
      }
      E+=p*Eik(tks[i],*k);
      sump+=p;
    }else{
        cout << "        ";
    }
      }
      cout << endl;
    }
    cout << endl << "T=" << 1/beta  << " E=" << E << " n="<< y.size() << "  F= " << F <<  endl <<  "----------" << endl;
  }
}

double DAClusterizerInZ::Eik	(	const track_t &	t,
		const vertex_t &	k
	)		const

Definition at line 63 of file DAClusterizerInZ.cc.

References DAClusterizerInZ::track_t::dz2, funct::pow(), DAClusterizerInZ::track_t::z, and DAClusterizerInZ::vertex_t::z.

                                                                      {
  return pow(t.z-k.z,2)/t.dz2;
}

vector< DAClusterizerInZ::track_t > DAClusterizerInZ::fill

(

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

tracks

)

const

Definition at line 20 of file DAClusterizerInZ.cc.

References DAClusterizerInZ::track_t::dz2, Measurement1D::error(), funct::exp(), DAClusterizerInZ::track_t::pi, position, funct::pow(), matplotRender::t, funct::tan(), DAClusterizerInZ::track_t::tt, Measurement1D::value(), DAClusterizerInZ::track_t::z, and DAClusterizerInZ::track_t::Z.

                    {
  // prepare track data for clustering
  vector<track_t> tks;
  for(vector<reco::TransientTrack>::const_iterator it=tracks.begin(); it!=tracks.end(); it++){
    track_t t;
    t.z=((*it).stateAtBeamLine().trackStateAtPCA()).position().z();
    double tantheta=tan(((*it).stateAtBeamLine().trackStateAtPCA()).momentum().theta());
    //  get the beam-spot
    reco::BeamSpot beamspot=(it->stateAtBeamLine()).beamSpot();
    t.dz2= pow((*it).track().dzError(),2)          // track errror
      + (pow(beamspot.BeamWidthX(),2)+pow(beamspot.BeamWidthY(),2))/pow(tantheta,2)  // beam-width induced
      + pow(vertexSize_,2);                        // intrinsic vertex size, safer for outliers and short lived decays
    if (d0CutOff_>0){
      Measurement1D IP=(*it).stateAtBeamLine().transverseImpactParameter();// error constains beamspot
      t.pi=1./(1.+exp(pow(IP.value()/IP.error(),2)-pow(d0CutOff_ ,2)));  // reduce weight for high ip tracks  
    }else{
      t.pi=1.;
    }
    t.tt=&(*it);
    t.Z=1.;
    tks.push_back(t);
  }
  return tks;
}

bool DAClusterizerInZ::merge	(	std::vector< vertex_t > &	y,
		int	nt
	)		const

Definition at line 210 of file DAClusterizerInZ.cc.

References gen::k, and detailsBasic3DVector::z.

                                                             {
  // merge clusters that collapsed or never separated, return true if vertices were merged, false otherwise

 if(y.size()<2)  return false;

  for(vector<vertex_t>::iterator k=y.begin(); (k+1)!=y.end(); k++){
    //if ((k+1)->z - k->z<1.e-2){  // note, no fabs here, maintains z-ordering  (with split()+merge() at every temperature)
    if (fabs((k+1)->z - k->z)<1.e-2){  // with fabs if only called after freeze-out (splitAll() at highter T)
      k->pk += (k+1)->pk;
      k->z=0.5*(k->z+(k+1)->z);
      y.erase(k+1);
      return true;  
    }
  }
  
  return false;
}

bool DAClusterizerInZ::purge	(	std::vector< vertex_t > &	y,
		std::vector< track_t > &	tks,
		double &	rho0,
		const double	beta
	)		const

Definition at line 230 of file DAClusterizerInZ.cc.

References gather_cfg::cout, funct::exp(), i, gen::k, reco::ParticleMasses::k0, nt, L1TEmulatorMonitor_cff::p, pi, and Gflash::Z.

                                                                                                              {
  // eliminate clusters with only one significant/unique track
  if(y.size()<2)  return false;
  
  unsigned int nt=tks.size();
  double sumpmin=nt;
  vector<vertex_t>::iterator k0=y.end();
  for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){ 
    int nUnique=0;
    double sump=0;
    double pmax=k->pk/(k->pk+rho0*exp(-beta*dzCutOff_*dzCutOff_));
    for(unsigned int i=0; i<nt; i++){
      if(tks[i].Z>0){
    //double p=pik(beta,tks[i],*k);
    double p=k->pk * exp(-beta*Eik(tks[i],*k)) / tks[i].Z;
    sump+=p;
    if( (p > 0.9*pmax) && (tks[i].pi>0) ){ nUnique++; }
      }
    }

    if((nUnique<2)&&(sump<sumpmin)){
      sumpmin=sump;
      k0=k;
    }
  }
 
  if(k0!=y.end()){
    if(verbose_){cout << "eliminating prototype at " << k0->z << " with sump=" << sumpmin << endl;}
    //rho0+=k0->pk;
    y.erase(k0);
    return true;
  }else{
    return false;
  }
}

void DAClusterizerInZ::splitAll	(	std::vector< track_t > &	tks,
		std::vector< vertex_t > &	y
	)		const

Definition at line 314 of file DAClusterizerInZ.cc.

References epsilon, gen::k, DAClusterizerInZ::vertex_t::pk, and DAClusterizerInZ::vertex_t::z.

                   {


  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)
  vector<vertex_t> y1;

  for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
    if ( ( (k==y.begin())|| (k-1)->z < k->z - zsep) && (((k+1)==y.end()  )|| (k+1)->z > k->z + zsep)) { 
      // isolated prototype, split
      vertex_t vnew;
      vnew.z  = k->z -epsilon;
      (*k).z  = k->z+epsilon;
      vnew.pk= 0.5* (*k).pk;
      (*k).pk= 0.5* (*k).pk;
      y1.push_back(vnew);
      y1.push_back(*k);

    }else if( y1.empty() || (y1.back().z < k->z -zsep)){
      y1.push_back(*k);
    }else{
      y1.back().z -=epsilon;
      k->z+=epsilon;
      y1.push_back(*k);
    }
  }// vertex loop
  
  y=y1;
}

double DAClusterizerInZ::update	(	double	beta,
		std::vector< track_t > &	tks,
		std::vector< vertex_t > &	y
	)		const

Definition at line 70 of file DAClusterizerInZ.cc.

References gather_cfg::cout, delta, funct::exp(), i, gen::k, nt, funct::pow(), and Gflash::Z.

Referenced by python.Vispa.Gui.VispaWidget.VispaWidget::autosize(), python.Vispa.Views.LineDecayView.LineDecayContainer::createObject(), python.Vispa.Views.LineDecayView.LineDecayContainer::deselectAllObjects(), python.Vispa.Gui.VispaWidgetOwner.VispaWidgetOwner::deselectAllWidgets(), python.Vispa.Gui.VispaWidget.VispaWidget::enableAutosizing(), python.Vispa.Gui.MenuWidget.MenuWidget::leaveEvent(), python.Vispa.Gui.VispaWidgetOwner.VispaWidgetOwner::mouseMoveEvent(), python.Vispa.Gui.MenuWidget.MenuWidget::mouseMoveEvent(), python.Vispa.Views.LineDecayView.LineDecayContainer::mouseMoveEvent(), python.Vispa.Gui.VispaWidgetOwner.VispaWidgetOwner::mouseReleaseEvent(), python.Vispa.Views.LineDecayView.LineDecayContainer::objectMoved(), python.Vispa.Views.LineDecayView.LineDecayContainer::removeObject(), python.Vispa.Gui.ConnectableWidget.ConnectableWidget::removePorts(), python.Vispa.Gui.FindDialog.FindDialog::reset(), python.Vispa.Gui.PortConnection.PointToPointConnection::select(), python.Vispa.Gui.VispaWidget.VispaWidget::select(), python.Vispa.Views.LineDecayView.LineDecayContainer::select(), python.Vispa.Gui.VispaWidget.VispaWidget::setText(), python.Vispa.Gui.VispaWidget.VispaWidget::setTitle(), python.Vispa.Gui.ZoomableWidget.ZoomableWidget::setZoom(), python.Vispa.Views.LineDecayView.LineDecayContainer::setZoom(), and python.Vispa.Gui.PortConnection.PointToPointConnection::updateConnection().

                             {
  //update weights and vertex positions
  // mass constrained annealing without noise
  // returns the squared sum of changes of vertex positions

  //cout << endl << endl << "DAClusterizerInZ::update" << endl;
  unsigned int nt=tks.size();

  //initialize sums
  double sumpi=0;
  for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
    k->se=0;    k->sw=0;   k->swz=0;
  }


  // loop over tracks
  for(unsigned int i=0; i<nt; i++){

    // update pik and Zi
    double Zi=0.;
    for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
      k->ei=exp(-beta*Eik(tks[i],*k));// cache exponential for one track at a time
      Zi += k->pk * k->ei;
    }
    tks[i].Z=Zi;

    // normalization for pk
    if (tks[i].Z>0){
      sumpi += tks[i].pi;
      // accumulate weighted z and weights for vertex update
      for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
    k->se += tks[i].pi* k->ei / Zi;
    double w = k->pk * tks[i].pi* k->ei / Zi / tks[i].dz2;
    k->sw  += w;
    k->swz += w * tks[i].z;
      }
    }else{
      sumpi += tks[i].pi;
    }
    

  } // end of track loop


  // now update z and pk
  double delta=0;
  for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
    if ( k->sw > 0){
      double znew=k->swz/k->sw; 
      delta+=pow(k->z-znew,2);
      k->z=znew;
    }else{
      edm::LogInfo("sumw") <<  "invalid sum of weights in fit: " << k->sw << endl;
      if(verbose_){cout << " a cluster melted away ?  pk=" << k->pk <<  " sumw=" << k->sw <<  endl;}
    }

    k->pk = k->pk * k->se / sumpi;
  }

  // return how much the prototypes moved
  return delta;
}

double DAClusterizerInZ::update	(	double	beta,
		std::vector< track_t > &	tks,
		std::vector< vertex_t > &	y,
		double &	rho0
	)		const

Definition at line 141 of file DAClusterizerInZ.cc.

References gather_cfg::cout, delta, funct::exp(), i, gen::k, nt, funct::pow(), and Gflash::Z.

Referenced by python.Vispa.Gui.VispaWidget.VispaWidget::autosize(), python.Vispa.Views.LineDecayView.LineDecayContainer::createObject(), python.Vispa.Views.LineDecayView.LineDecayContainer::deselectAllObjects(), python.Vispa.Gui.VispaWidgetOwner.VispaWidgetOwner::deselectAllWidgets(), python.Vispa.Gui.VispaWidget.VispaWidget::enableAutosizing(), python.Vispa.Gui.MenuWidget.MenuWidget::leaveEvent(), python.Vispa.Gui.VispaWidgetOwner.VispaWidgetOwner::mouseMoveEvent(), python.Vispa.Gui.MenuWidget.MenuWidget::mouseMoveEvent(), python.Vispa.Views.LineDecayView.LineDecayContainer::mouseMoveEvent(), python.Vispa.Gui.VispaWidgetOwner.VispaWidgetOwner::mouseReleaseEvent(), python.Vispa.Views.LineDecayView.LineDecayContainer::objectMoved(), python.Vispa.Views.LineDecayView.LineDecayContainer::removeObject(), python.Vispa.Gui.ConnectableWidget.ConnectableWidget::removePorts(), python.Vispa.Gui.FindDialog.FindDialog::reset(), python.Vispa.Gui.PortConnection.PointToPointConnection::select(), python.Vispa.Gui.VispaWidget.VispaWidget::select(), python.Vispa.Views.LineDecayView.LineDecayContainer::select(), python.Vispa.Gui.VispaWidget.VispaWidget::setText(), python.Vispa.Gui.VispaWidget.VispaWidget::setTitle(), python.Vispa.Gui.ZoomableWidget.ZoomableWidget::setZoom(), python.Vispa.Views.LineDecayView.LineDecayContainer::setZoom(), and python.Vispa.Gui.PortConnection.PointToPointConnection::updateConnection().

                      {
  // MVF style, no more vertex weights, update tracks weights and vertex positions, with noise 
  // returns the squared sum of changes of vertex positions

  unsigned int nt=tks.size();

  //initialize sums
  for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
    k->se=0;    k->sw=0;   k->swz=0;
  }


  // loop over tracks
  for(unsigned int i=0; i<nt; i++){

    // update pik and Zi
    double Zi=rho0*exp(-beta*dzCutOff_*dzCutOff_);// cut-off
    for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
      k->ei=exp(-beta*Eik(tks[i],*k));// cache exponential for one track at a time
      Zi += k->pk * k->ei;
    }
    tks[i].Z=Zi;

    // normalization
    if (tks[i].Z>0){
      //double pi0=exp(-beta*dzCutOff_*dzCutOff_)/tks[i].Z;
      //sumpi += tks[i].pi*(1.-pi0*rho0);  // exclude rho0 from the normalization of the cluster "mass"
       // accumulate weighted z and weights for vertex update
      for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
    k->se += tks[i].pi* k->ei / Zi;
    double w = k->pk * tks[i].pi* k->ei / Zi / tks[i].dz2;
    k->sw  += w;
    k->swz += w * tks[i].z;
      }
    }
    //else{  sumpi += tks[i].pi;  }


  } // end of track loop


  // now update z and pk
  double delta=0;
  for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
    if ( k->sw > 0){
      double znew=k->swz/k->sw; 
      delta+=pow(k->z-znew,2);
      k->z=znew;
    }else{
      edm::LogInfo("sumw") <<  "invalid sum of weights in fit: " << k->sw << endl;
      if(verbose_){cout << " a cluster melted away ?  pk=" << k->pk <<  " sumw=" << k->sw <<  endl;}
    }

    //k->pk *= k->se / sumpi; 
  }

  // return how much the prototypes moved
  return delta;
}

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

Definition at line 449 of file DAClusterizerInZ.cc.

References beta, gather_cfg::cout, hcal_timing_source_file_cfg::dump, funct::exp(), i, gen::k, relval_steps::merge(), nt, L1TEmulatorMonitor_cff::p, pi, DAClusterizerInZ::vertex_t::pk, pos, update, v, detailsBasic3DVector::y, DAClusterizerInZ::vertex_t::z, and Gflash::Z.

{
 
  vector<track_t> tks=fill(tracks);
  unsigned int nt=tracks.size();
  double rho0=0.0;  // start with no outlier rejection

  vector< TransientVertex > clusters;
  if (tks.empty()) return clusters;

  vector<vertex_t> y; // the vertex prototypes

  // initialize:single vertex at infinite temperature
  vertex_t vstart;
  vstart.z=0.;
  vstart.pk=1.;
  y.push_back(vstart);
  int niter=0;      // number of iterations
  

  // estimate first critical temperature
  double beta=beta0(betamax_, tks, y);
  niter=0; while((update(beta, tks,y)>1.e-6)  && (niter++ < maxIterations_)){ }

  
 // annealing loop, stop when T<Tmin  (i.e. beta>1/Tmin)
  while(beta<betamax_){ 

    beta=beta/coolingFactor_;
    splitAll(tks,y);

    // make sure we are not too far from equilibrium before cooling further
    niter=0; while((update(beta, tks,y)>1.e-6)  && (niter++ < maxIterations_)){ }

  }


  // merge collapsed clusters 
  while(merge(y,tks.size())){} 
  if(verbose_  ){ cout << "dump after 1st merging " << endl;  dump(beta,y,tks,2);}



  // switch on outlier rejection
  //rho0=exp(beta*dzCutOff_*dzCutOff_)/nt; if(rho0>0.1) rho0=0.1;
  rho0=1./nt; for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){ k->pk =1.; }  // democratic
  niter=0; while((update(beta, tks,y,rho0) > 1.e-8)  && (niter++ < maxIterations_)){  }
  if(verbose_  ){ cout << "rho0=" << rho0 <<   " niter=" << niter <<  endl; dump(beta,y,tks,2);}

  
  // merge again  (some cluster split by outliers collapse here)
  while(merge(y,tks.size())){}  
  if(verbose_  ){ cout << "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,rho0) > 1.e-6)  && (niter++ < maxIterations_)){  }
    } 
    beta/=coolingFactor_;
    niter=0; while((update(beta, tks,y,rho0) > 1.e-6)  && (niter++ < maxIterations_)){  }
  }

  if(verbose_){
   cout << "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
  for(unsigned int i=0; i<nt; i++){  
    tks[i].Z=rho0*exp(-beta*dzCutOff_*dzCutOff_);
    for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){ 
      tks[i].Z+=k->pk * exp(-beta*Eik(tks[i],*k));
    }
  }


  for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){ 
    GlobalPoint pos(0, 0, k->z);
    vector< reco::TransientTrack > vertexTracks;
    for(unsigned int i=0; i<nt; i++){
      if(tks[i].Z>0){
    double p=k->pk * exp(-beta*Eik(tks[i],*k)) / tks[i].Z;
    if( (tks[i].pi>0) && ( p > 0.5 ) ){ vertexTracks.push_back(*(tks[i].tt)); tks[i].Z=0; } // setting Z=0 excludes double assignment
      }
    }
    TransientVertex v(pos, dummyError, vertexTracks, 0);
    clusters.push_back(v);
  }


  return clusters;

}

Member Data Documentation

float DAClusterizerInZ::betamax_

private

Definition at line 96 of file DAClusterizerInZ.h.

float DAClusterizerInZ::betastop_

private

Definition at line 97 of file DAClusterizerInZ.h.

double DAClusterizerInZ::coolingFactor_

private

Definition at line 95 of file DAClusterizerInZ.h.

double DAClusterizerInZ::d0CutOff_

private

Definition at line 99 of file DAClusterizerInZ.h.

double DAClusterizerInZ::dzCutOff_

private

Definition at line 98 of file DAClusterizerInZ.h.

int DAClusterizerInZ::maxIterations_

private

Definition at line 94 of file DAClusterizerInZ.h.

bool DAClusterizerInZ::verbose_

private

Definition at line 92 of file DAClusterizerInZ.h.

float DAClusterizerInZ::vertexSize_

private

Definition at line 93 of file DAClusterizerInZ.h.