CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZ_vect.cc

Go to the documentation of this file.
00001 #include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ_vect.h"
00002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00003 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
00004 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
00005 
00006 #include <cmath>
00007 #include <cassert>
00008 #include <limits>
00009 #include <iomanip>
00010 
00011 #include "vdt/vdtMath.h"
00012 
00013 using namespace std;
00014 
00015 DAClusterizerInZ_vect::DAClusterizerInZ_vect(const edm::ParameterSet& conf) {
00016 
00017         // some defaults to avoid uninitialized variables
00018         verbose_ = conf.getUntrackedParameter<bool> ("verbose", false);
00019         use_vdt_ = conf.getUntrackedParameter<bool> ("use_vdt", false);
00020 
00021         betamax_ = 0.1;
00022         betastop_ = 1.0;
00023         coolingFactor_ = 0.6;
00024         maxIterations_ = 100;
00025         vertexSize_ = 0.01; // 0.1 mm
00026         dzCutOff_ = 4.0;  
00027 
00028         // configure
00029 
00030         double Tmin = conf.getParameter<double> ("Tmin");
00031         vertexSize_ = conf.getParameter<double> ("vertexSize");
00032         coolingFactor_ = conf.getParameter<double> ("coolingFactor");
00033         useTc_=true;
00034         if(coolingFactor_<0){
00035           coolingFactor_=-coolingFactor_; 
00036           useTc_=false;
00037         }
00038         d0CutOff_ = conf.getParameter<double> ("d0CutOff");
00039         dzCutOff_ = conf.getParameter<double> ("dzCutOff");
00040         maxIterations_ = 100;
00041         if (Tmin == 0) {
00042                 LogDebug("DAClusterizerinZ_vectorized")  << "DAClusterizerInZ: invalid Tmin" << Tmin
00043                                 << "  reset do default " << 1. / betamax_ << endl;
00044         } else {
00045                 betamax_ = 1. / Tmin;
00046         }
00047 
00048 }
00049 
00050 inline double DAClusterizerInZ_vect::local_exp( double const& inp) const
00051                 {
00052                         if ( use_vdt_)
00053                           return vdt::fast_exp( inp );
00054 
00055                         return std::exp( inp );
00056                 }
00057 
00058 inline void DAClusterizerInZ_vect::local_exp_list( double* arg_inp, double* arg_out,const int arg_arr_size) const
00059                 {
00060                         if ( use_vdt_){
00061 //                          std::cout << "I use vdt!\n";
00062                             for(int i=0; i!=arg_arr_size; ++i ) arg_out[i]=vdt::fast_exp(arg_inp[i]);
00063                            // vdt::fast_expv(arg_arr_size, arg_inp, arg_out);
00064                         }
00065                         else
00066                           for(int i=0; i!=arg_arr_size; ++i ) arg_out[i]=std::exp(arg_inp[i]);
00067                 }
00068 
00069 //todo: use r-value possibility of c++11 here
00070 DAClusterizerInZ_vect::track_t DAClusterizerInZ_vect::fill(const vector<
00071                 reco::TransientTrack> & tracks) const {
00072 
00073         // prepare track data for clustering
00074         track_t tks;
00075         for (vector<reco::TransientTrack>::const_iterator it = tracks.begin(); it
00076                         != tracks.end(); it++)
00077         {
00078                 double t_pi;
00079                 double t_z = ((*it).stateAtBeamLine().trackStateAtPCA()).position().z();
00080                 double phi=((*it).stateAtBeamLine().trackStateAtPCA()).momentum().phi();
00081                 double tantheta = tan(
00082                                 ((*it).stateAtBeamLine().trackStateAtPCA()).momentum().theta());
00083                 //  get the beam-spot
00084                 reco::BeamSpot beamspot = (it->stateAtBeamLine()).beamSpot();
00085                 double t_dz2 = 
00086                       pow((*it).track().dzError(), 2) // track errror
00087                   + (pow(beamspot.BeamWidthX()*cos(phi),2)+pow(beamspot.BeamWidthY()*sin(phi),2))/pow(tantheta,2)  // beam-width
00088                   + pow(vertexSize_, 2); // intrinsic vertex size, safer for outliers and short lived decays
00089                 if (d0CutOff_ > 0) {
00090                         Measurement1D IP =
00091                                         (*it).stateAtBeamLine().transverseImpactParameter();// error constains beamspot
00092                         t_pi = 1. / (1. + local_exp(pow(IP.value() / IP.error(), 2) - pow(
00093                                         d0CutOff_, 2))); // reduce weight for high ip tracks
00094                 } else {
00095                         t_pi = 1.;
00096                 }
00097 
00098                 tks.AddItem(t_z, t_dz2, &(*it), t_pi);
00099         }
00100         tks.ExtractRaw();
00101 
00102         if (verbose_) {
00103                 LogDebug("DAClusterizerinZ_vectorized") << "Track count " << tks.GetSize() << std::endl;
00104         }
00105 
00106         return tks;
00107 }
00108 
00109 
00110 double DAClusterizerInZ_vect::Eik(double const& t_z, double const& k_z, double const& t_dz2) const
00111 {
00112         return pow(t_z - k_z, 2) / t_dz2;
00113 }
00114 
00115 
00116 double DAClusterizerInZ_vect::update(double beta, track_t & gtracks,
00117                 vertex_t & gvertices, bool useRho0, double & rho0) const {
00118 
00119         //update weights and vertex positions
00120         // mass constrained annealing without noise
00121         // returns the squared sum of changes of vertex positions
00122 
00123         const unsigned int nt = gtracks.GetSize();
00124         const unsigned int nv = gvertices.GetSize();
00125 
00126         //initialize sums
00127         double sumpi = 0;
00128 
00129         // to return how much the prototype moved
00130         double delta = 0;
00131 
00132         unsigned int itrack;
00133         unsigned int ivertex;
00134 
00135         // intial value of a sum
00136         double Z_init = 0;
00137 
00138         // define kernels
00139         auto kernel_calc_exp_arg = [ &beta, nv ] (
00140                         const unsigned int itrack,
00141                         track_t const& tracks,
00142                         vertex_t const& vertices )
00143         {
00144                 const double track_z = tracks._z[itrack];
00145                 const double track_dz2 = tracks._dz2[itrack];
00146 
00147                 // auto-vectorized
00148                 for ( unsigned int ivertex = 0; ivertex < nv; ++ivertex)
00149                 {
00150                         double mult_res = ( track_z - vertices._z[ivertex] );
00151                         vertices._ei_cache[ivertex] = -beta *
00152                         ( mult_res * mult_res ) / track_dz2;
00153                 }
00154         };
00155 
00156         // auto kernel_add_Z = [ nv, &Z_init ] (vertex_t const& vertices) raises a warning
00157         // we declare the return type with " -> double"
00158         auto kernel_add_Z = [ nv, &Z_init ] (vertex_t const& vertices) -> double
00159         {
00160                 double ZTemp = Z_init;
00161  
00162                 for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
00163 
00164                         ZTemp += vertices._pk[ivertex] * vertices._ei[ivertex];
00165                 }
00166                 return ZTemp;
00167         };
00168 
00169         auto kernel_calc_normalization = [ &beta, nv ] (const unsigned int track_num,
00170                         track_t & tks_vec,
00171                         vertex_t & y_vec )
00172         {
00173                 double w;
00174 
00175                 const double tmp_trk_pi = tks_vec._pi[track_num];
00176                 const double tmp_trk_Z_sum = tks_vec._Z_sum[track_num];
00177                 const double tmp_trk_dz2 = tks_vec._dz2[track_num];
00178                 const double tmp_trk_z = tks_vec._z[track_num];
00179 
00180                 // auto-vectorized
00181                 for (unsigned int k = 0; k < nv; ++k) {
00182                         y_vec._se[k] += tmp_trk_pi* y_vec._ei[k] / tmp_trk_Z_sum;
00183                         w = y_vec._pk[k] * tmp_trk_pi * y_vec._ei[k] / tmp_trk_Z_sum / tmp_trk_dz2;
00184                         y_vec._sw[k]  += w;
00185                         y_vec._swz[k] += w * tmp_trk_z;
00186                         y_vec._swE[k] += w * y_vec._ei_cache[k]/(-beta);
00187                 }
00188         };
00189 
00190         // not vectorized
00191         for (ivertex = 0; ivertex < nv; ++ivertex) {
00192                 gvertices._se[ivertex] = 0.0;
00193                 gvertices._sw[ivertex] = 0.0;
00194                 gvertices._swz[ivertex] = 0.0;
00195                 gvertices._swE[ivertex] = 0.0;
00196         }
00197 
00198 
00199         // independpent of loop
00200         if ( useRho0 )
00201         {
00202                 Z_init = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_); // cut-off
00203         }
00204 
00205         // loop over tracks
00206         for (itrack = 0; itrack < nt; ++itrack) {
00207                 kernel_calc_exp_arg(itrack, gtracks, gvertices);
00208                 local_exp_list(gvertices._ei_cache, gvertices._ei, nv);
00209                 //vdt::cephes_single_exp_vect( y_vec._ei, nv );
00210 
00211                 gtracks._Z_sum[itrack] = kernel_add_Z(gvertices);
00212 
00213                 // used in the next major loop to follow
00214                 if (!useRho0)
00215                         sumpi += gtracks._pi[itrack];
00216 
00217                 if (gtracks._Z_sum[itrack] > 0) {
00218                         kernel_calc_normalization(itrack, gtracks, gvertices);
00219                 }
00220         }
00221 
00222         // now update z and pk
00223         auto kernel_calc_z = [ &delta, &sumpi, nv, this, useRho0 ] (vertex_t & vertices )
00224         {
00225                 // does not vectorizes
00226                 for (unsigned int ivertex = 0; ivertex < nv; ++ ivertex )
00227                 {
00228                         if (vertices._sw[ivertex] > 0)
00229                         {
00230                                 double znew = vertices._swz[ ivertex ] / vertices._sw[ ivertex ];
00231 
00232                                 // prevents from vectorizing
00233                                 delta += pow( vertices._z[ ivertex ] - znew, 2 );
00234                                 vertices._z[ ivertex ] = znew;
00235                         }
00236                         else {
00237                                 edm::LogInfo("sumw") << "invalid sum of weights in fit: " << vertices._sw[ivertex]
00238                                 << endl;
00239                                 if (this->verbose_) {
00240                                         LogDebug("DAClusterizerinZ_vectorized")  << " a cluster melted away ?  pk=" << vertices._pk[ ivertex ] << " sumw="
00241                                         << vertices._sw[ivertex] << endl;
00242                                 }
00243                         }
00244 
00245                         // dont do, if rho cut
00246                         if ( ! useRho0 )
00247                         {
00248                                 vertices._pk[ ivertex ] = vertices._pk[ ivertex ] * vertices._se[ ivertex ] / sumpi;
00249                         }
00250                 }
00251         };
00252 
00253         kernel_calc_z(gvertices);
00254 
00255         // return how much the prototypes moved
00256         return delta;
00257 }
00258 
00259 
00260 
00261 bool DAClusterizerInZ_vect::merge(vertex_t & y, double & beta)const{
00262   // merge clusters that collapsed or never separated,
00263   // only merge if the estimated critical temperature of the merged vertex is below the current temperature
00264   // return true if vertices were merged, false otherwise
00265   const unsigned int nv = y.GetSize();
00266 
00267   if (nv < 2)
00268     return false;
00269 
00270 
00271   for (unsigned int k = 0; (k + 1) < nv; k++) {
00272     if (fabs(y._z[k + 1] - y._z[k]) < 2.e-2) {
00273       double rho=y._pk[k] + y._pk[k+1];
00274       double swE=y._swE[k]+y._swE[k+1] - y._pk[k]*y._pk[k+1] /rho *pow(y._z[k+1]-y._z[k],2);
00275       double Tc=2*swE/(y._sw[k]+y._sw[k+1]);
00276 
00277       if(Tc*beta<1){
00278         if(rho>0){
00279           y._z[k] = (y._pk[k]*y._z[k] + y._pk[k+1]*y._z[k + 1])/rho;
00280         }else{
00281           y._z[k] = 0.5 * (y._z[k] + y._z[k + 1]);
00282         }
00283         y._pk[k] = rho;
00284         y._sw[k]+=y._sw[k+1];
00285         y._swE[k]=swE;
00286         y.RemoveItem(k+1);
00287         return true;
00288       }
00289     }
00290   }
00291 
00292   return false;
00293 }
00294 
00295 
00296 
00297 
00298 
00299 bool DAClusterizerInZ_vect::merge(vertex_t & y) const {
00300         // merge clusters that collapsed or never separated, return true if vertices were merged, false otherwise
00301 
00302         const unsigned int nv = y.GetSize();
00303 
00304         if (nv < 2)
00305                 return false;
00306 
00307         for (unsigned int k = 0; (k + 1) < nv; k++) {
00308                 //if ((k+1)->z - k->z<1.e-2){  // note, no fabs here, maintains z-ordering  (with split()+merge() at every temperature)
00309                 if (fabs(y._z[k + 1] - y._z[k]) < 1.e-2) { // with fabs if only called after freeze-out (splitAll() at highter T)
00310                         y._pk[k] += y._pk[k + 1];
00311                         y._z[k] = 0.5 * (y._z[k] + y._z[k + 1]);
00312 
00313                         y.RemoveItem(k + 1);
00314 
00315                         if ( verbose_)
00316                         {
00317                                 LogDebug("DAClusterizerinZ_vectorized") << "Merging vertices k = " << k << std::endl;
00318                         }
00319 
00320                         return true;
00321                 }
00322         }
00323 
00324         return false;
00325 }
00326 
00327 bool DAClusterizerInZ_vect::purge(vertex_t & y, track_t & tks, double & rho0,
00328                 const double beta) const {
00329         // eliminate clusters with only one significant/unique track
00330         const unsigned int nv = y.GetSize();
00331         const unsigned int nt = tks.GetSize();
00332 
00333         if (nv < 2)
00334                 return false;
00335 
00336         double sumpmin = nt;
00337         unsigned int k0 = nv;
00338 
00339         for (unsigned int k = 0; k < nv; k++) {
00340 
00341                 int nUnique = 0;
00342                 double sump = 0;
00343                 double pmax = y._pk[k] / (y._pk[k] + rho0 * local_exp(-beta * dzCutOff_
00344                                 * dzCutOff_));
00345 
00346                 for (unsigned int i = 0; i < nt; i++) {
00347 
00348                         if (tks._Z_sum[i] > 0) {
00349                                 //double p=pik(beta,tks[i],*k);
00350                                 double p = y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k],
00351                                                 tks._dz2[i])) / tks._Z_sum[i];
00352                                 sump += p;
00353                                 if ((p > 0.9 * pmax) && (tks._pi[i] > 0)) {
00354                                         nUnique++;
00355                                 }
00356                         }
00357                 }
00358 
00359                 if ((nUnique < 2) && (sump < sumpmin)) {
00360                         sumpmin = sump;
00361                         k0 = k;
00362                 }
00363         }
00364 
00365         if (k0 != nv) {
00366                 if (verbose_) {
00367                         LogDebug("DAClusterizerinZ_vectorized")  << "eliminating prototype at " << y._z[k0] << " with sump="
00368                                         << sumpmin << endl;
00369                 }
00370                 //rho0+=k0->pk;
00371                 y.RemoveItem(k0);
00372                 return true;
00373         } else {
00374                 return false;
00375         }
00376 }
00377 
00378 double DAClusterizerInZ_vect::beta0(double betamax, track_t & tks, vertex_t & y) const {
00379 
00380         double T0 = 0; // max Tc for beta=0
00381         // estimate critical temperature from beta=0 (T=inf)
00382         const unsigned int nt = tks.GetSize();
00383         const unsigned int nv = y.GetSize();
00384 
00385         for (unsigned int k = 0; k < nv; k++) {
00386 
00387                 // vertex fit at T=inf
00388                 double sumwz = 0;
00389                 double sumw = 0;
00390                 for (unsigned int i = 0; i < nt; i++) {
00391                         double w = tks._pi[i] / tks._dz2[i];
00392                         sumwz += w * tks._z[i];
00393                         sumw += w;
00394                 }
00395                 y._z[k] = sumwz / sumw;
00396 
00397                 // estimate Tcrit, eventually do this in the same loop
00398                 double a = 0, b = 0;
00399                 for (unsigned int i = 0; i < nt; i++) {
00400                         double dx = tks._z[i] - (y._z[k]);
00401                         double w = tks._pi[i] / tks._dz2[i];
00402                         a += w * pow(dx, 2) / tks._dz2[i];
00403                         b += w;
00404                 }
00405                 double Tc = 2. * a / b; // the critical temperature of this vertex
00406                 if (Tc > T0)
00407                         T0 = Tc;
00408         }// vertex loop (normally there should be only one vertex at beta=0)
00409 
00410         if (T0 > 1. / betamax) {
00411                 return betamax / pow(coolingFactor_, int(log(T0 * betamax) / log(
00412                                 coolingFactor_)) - 1);
00413         } else {
00414                 // ensure at least one annealing step
00415                 return betamax / coolingFactor_;
00416         }
00417 }
00418 
00419 
00420 bool DAClusterizerInZ_vect::split(const double beta,  track_t &tks, vertex_t & y ) const{
00421   // split only critical vertices (Tc >~ T=1/beta   <==>   beta*Tc>~1)
00422   // an update must have been made just before doing this (same beta, no merging)
00423   // returns true if at least one cluster was split
00424 
00425   double epsilon=1e-3;      // split all single vertices by 10 um
00426   unsigned int nv = y.GetSize();
00427 
00428   // avoid left-right biases by splitting highest Tc first
00429 
00430   std::vector<std::pair<double, unsigned int> > critical;
00431   for(unsigned int k=0; k<nv; k++){
00432     double Tc= 2*y._swE[k]/y._sw[k];
00433     if (beta*Tc > 1.){
00434       critical.push_back( make_pair(Tc, k));
00435     }
00436   }
00437   if (critical.size()==0) return false;
00438   stable_sort(critical.begin(), critical.end(), std::greater<std::pair<double, unsigned int> >() );
00439 
00440 
00441   bool split=false;
00442   const unsigned int nt = tks.GetSize();
00443 
00444   for(unsigned int ic=0; ic<critical.size(); ic++){
00445     unsigned int k=critical[ic].second;
00446     // estimate subcluster positions and weight
00447     double p1=0, z1=0, w1=0;
00448     double p2=0, z2=0, w2=0;
00449     for(unsigned int i=0; i<nt; i++){
00450       if (tks._Z_sum[i] > 0) {
00451         double p = y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k],
00452                                                     tks._dz2[i])) / tks._Z_sum[i];
00453 
00454         double w=p/tks._dz2[i];
00455         if(tks._z[i]<y._z[k]){
00456           p1+=p; z1+=w*tks._z[i]; w1+=w;
00457         }else{
00458           p2+=p; z2+=w*tks._z[i]; w2+=w;
00459         }
00460       }
00461     }
00462     if(w1>0){  z1=z1/w1;} else{z1=y._z[k]-epsilon;}
00463     if(w2>0){  z2=z2/w2;} else{z2=y._z[k]+epsilon;}
00464 
00465     // reduce split size if there is not enough room
00466     if( ( k   > 0 ) && ( y._z[k-1]>=z1 ) ){ z1=0.5*(y._z[k]+y._z[k-1]); }
00467     if( ( k+1 < nv) && ( y._z[k+1]<=z2 ) ){ z2=0.5*(y._z[k]+y._z[k+1]); }
00468 
00469     // split if the new subclusters are significantly separated
00470     if( (z2-z1)>epsilon){
00471       split=true;
00472       double pk1=p1*y._pk[k]/(p1+p2);
00473       double pk2=p2*y._pk[k]/(p1+p2);
00474       y._z[k]  =  z2;
00475       y._pk[k] = pk2;
00476       y.InsertItem(k, z1, pk1);
00477       nv++;
00478 
00479      // adjust remaining pointers
00480       for(unsigned int jc=ic; jc<critical.size(); jc++){
00481         if (critical[jc].second>k) {critical[jc].second++;}
00482       }
00483     }
00484   }
00485   return split;
00486 }
00487 
00488 
00489 
00490 void DAClusterizerInZ_vect::splitAll( vertex_t & y) const {
00491 
00492         const unsigned int nv = y.GetSize();
00493 
00494         double epsilon = 1e-3; // split all single vertices by 10 um
00495         double zsep = 2 * epsilon; // split vertices that are isolated by at least zsep (vertices that haven't collapsed)
00496         vertex_t y1;
00497 
00498         if (verbose_)
00499         {
00500                 LogDebug("DAClusterizerinZ_vectorized") << "Before Split "<< std::endl;
00501                 y.DebugOut();
00502         }
00503 
00504         for (unsigned int k = 0; k < nv; k++) {
00505 
00506                 if (
00507                                 ( (k == 0)              || ( y._z[k - 1]        < (y._z[k] - zsep)) ) &&
00508                         ( ((k + 1) == nv)       || ( y._z[k + 1]        > (y._z[k] + zsep)) )   )
00509                 {
00510                         // isolated prototype, split
00511                         double new_z = y.z[k] - epsilon;
00512                         y.z[k] = y.z[k] + epsilon;
00513 
00514                         double new_pk = 0.5 * y.pk[k];
00515                         y.pk[k] = 0.5 * y.pk[k];
00516 
00517                         y1.AddItem(new_z, new_pk);
00518                         y1.AddItem(y._z[k], y._pk[k]);
00519 
00520                 }
00521                 else if ( (y1.GetSize() == 0 ) ||
00522                                 (y1._z[y1.GetSize() - 1] <  (y._z[k] - zsep)  ))
00523                 {
00524 
00525                         y1.AddItem(y._z[k], y._pk[k]);
00526                 }
00527                 else
00528                 {
00529                         y1._z[y1.GetSize() - 1] = y1._z[y1.GetSize() - 1] - epsilon;
00530                         y._z[k] = y._z[k] + epsilon;
00531                         y1.AddItem( y._z[k] , y._pk[k]);
00532                 }
00533         }// vertex loop
00534 
00535         y = y1;
00536         y.ExtractRaw();
00537 
00538         if (verbose_)
00539         {
00540                 LogDebug("DAClusterizerinZ_vectorized") << "After split " << std::endl;
00541                 y.DebugOut();
00542         }
00543 }
00544 
00545 
00546 void DAClusterizerInZ_vect::dump(const double beta, const vertex_t & y,
00547                 const track_t & tks, int verbosity) const {
00548 
00549         const unsigned int nv = y.GetSize();
00550         const unsigned int nt = tks.GetSize();
00551 
00552         LogDebug("DAClusterizerinZ_vectorized")  << "-----DAClusterizerInZ::dump ----" << endl;
00553         LogDebug("DAClusterizerinZ_vectorized")  << "beta=" << beta << "   betamax= " << betamax_ << endl;
00554         LogDebug("DAClusterizerinZ_vectorized")  << "                                                               z= ";
00555         LogDebug("DAClusterizerinZ_vectorized")  << setprecision(4);
00556         for (unsigned int ivertex = 0; ivertex < nv; ++ ivertex) {
00557                 LogDebug("DAClusterizerinZ_vectorized")  << setw(8) << fixed << y._z[ivertex];
00558         }
00559         LogDebug("DAClusterizerinZ_vectorized")  << endl << "T=" << setw(15) << 1. / beta
00560                         << "                                            Tc= ";
00561         LogDebug("DAClusterizerinZ_vectorized")  << endl
00562                         << "                                                               pk=";
00563         double sumpk = 0;
00564         for (unsigned int ivertex = 0; ivertex < nv; ++ ivertex) {
00565                 LogDebug("DAClusterizerinZ_vectorized")  << setw(8) << setprecision(3) << fixed << y._pk[ivertex];
00566                 sumpk += y._pk[ivertex];
00567         }
00568         LogDebug("DAClusterizerinZ_vectorized")  << endl;
00569 
00570         if (verbosity > 0) {
00571                 double E = 0, F = 0;
00572                 LogDebug("DAClusterizerinZ_vectorized")  << endl;
00573                 LogDebug("DAClusterizerinZ_vectorized") 
00574                 << "----       z +/- dz                ip +/-dip       pt    phi  eta    weights  ----"
00575                 << endl;
00576                 LogDebug("DAClusterizerinZ_vectorized")  << setprecision(4);
00577                 for (unsigned int i = 0; i < nt; i++) {
00578                         if (tks._Z_sum[i] > 0) {
00579                                 F -= log(tks._Z_sum[i]) / beta;
00580                         }
00581                         double tz = tks._z[i];
00582                         LogDebug("DAClusterizerinZ_vectorized")  << setw(3) << i << ")" << setw(8) << fixed << setprecision(4)
00583                          << tz << " +/-" << setw(6) << sqrt(tks._dz2[i]);
00584 
00585                         if (tks.tt[i]->track().quality(reco::TrackBase::highPurity)) {
00586                                 LogDebug("DAClusterizerinZ_vectorized")  << " *";
00587                         } else {
00588                                 LogDebug("DAClusterizerinZ_vectorized")  << "  ";
00589                         }
00590                         if (tks.tt[i]->track().hitPattern().hasValidHitInFirstPixelBarrel()) {
00591                                 LogDebug("DAClusterizerinZ_vectorized")  << "+";
00592                         } else {
00593                                 LogDebug("DAClusterizerinZ_vectorized")  << "-";
00594                         }
00595                         LogDebug("DAClusterizerinZ_vectorized")  << setw(1)
00596                          << tks.tt[i]->track().hitPattern().pixelBarrelLayersWithMeasurement(); // see DataFormats/TrackReco/interface/HitPattern.h
00597                         LogDebug("DAClusterizerinZ_vectorized")  << setw(1)
00598                          << tks.tt[i]->track().hitPattern().pixelEndcapLayersWithMeasurement();
00599                         LogDebug("DAClusterizerinZ_vectorized")  << setw(1) << hex
00600                                         << tks.tt[i]->track().hitPattern().trackerLayersWithMeasurement()
00601                                         - tks.tt[i]->track().hitPattern().pixelLayersWithMeasurement()
00602                                         << dec;
00603                         LogDebug("DAClusterizerinZ_vectorized")  << "=" << setw(1) << hex
00604                                         << tks.tt[i]->track().trackerExpectedHitsOuter().numberOfHits()
00605                                         << dec;
00606 
00607                         Measurement1D IP =
00608                                         tks.tt[i]->stateAtBeamLine().transverseImpactParameter();
00609                         LogDebug("DAClusterizerinZ_vectorized")  << setw(8) << IP.value() << "+/-" << setw(6) << IP.error();
00610                         LogDebug("DAClusterizerinZ_vectorized")  << " " << setw(6) << setprecision(2)
00611                          << tks.tt[i]->track().pt() * tks.tt[i]->track().charge();
00612                         LogDebug("DAClusterizerinZ_vectorized")  << " " << setw(5) << setprecision(2)
00613                          << tks.tt[i]->track().phi() << " " << setw(5)
00614                          << setprecision(2) << tks.tt[i]->track().eta();
00615 
00616                         double sump = 0.;
00617                         for (unsigned int ivertex = 0; ivertex < nv; ++ ivertex) {
00618                                 if ((tks._pi[i] > 0) && (tks._Z_sum[i] > 0)) {
00619                                         //double p=pik(beta,tks[i],*k);
00620                                         double p = y._pk[ivertex]  * exp(-beta * Eik(tks._z[i], y._z[ivertex], tks._dz2[i])) / tks._Z_sum[i];
00621                                         if (p > 0.0001) {
00622                                                 LogDebug("DAClusterizerinZ_vectorized")  << setw(8) << setprecision(3) << p;
00623                                         } else {
00624                                                 LogDebug("DAClusterizerinZ_vectorized")  << "    .   ";
00625                                         }
00626                                         E += p * Eik(tks._z[i], y._z[ivertex], tks._dz2[i]);
00627                                         sump += p;
00628                                 } else {
00629                                         LogDebug("DAClusterizerinZ_vectorized")  << "        ";
00630                                 }
00631                         }
00632                         LogDebug("DAClusterizerinZ_vectorized")  << endl;
00633                 }
00634                 LogDebug("DAClusterizerinZ_vectorized")  << endl << "T=" << 1 / beta << " E=" << E << " n=" << y.GetSize()
00635                          << "  F= " << F << endl << "----------" << endl;
00636         }
00637 }
00638 
00639 vector<TransientVertex> DAClusterizerInZ_vect::vertices(const vector<
00640                 reco::TransientTrack> & tracks, const int verbosity) const {
00641 
00642 
00643         track_t tks = fill(tracks);
00644         tks.ExtractRaw();
00645 
00646         unsigned int nt = tracks.size();
00647         double rho0 = 0.0; // start with no outlier rejection
00648 
00649         vector<TransientVertex> clusters;
00650         if (tks.GetSize() == 0)
00651                 return clusters;
00652 
00653         vertex_t y; // the vertex prototypes
00654 
00655         // initialize:single vertex at infinite temperature
00656         y.AddItem( 0, 1.0);
00657 
00658         int niter = 0; // number of iterations
00659 
00660 
00661         // estimate first critical temperature
00662         double beta = beta0(betamax_, tks, y);
00663         if ( verbose_)
00664         {
00665                 LogDebug("DAClusterizerinZ_vectorized") << "Beta0 is " << beta << std::endl;
00666         }
00667 
00668         niter = 0;
00669         while ((update(beta, tks, y, false, rho0) > 1.e-6) && (niter++
00670                         < maxIterations_)) {
00671         }
00672 
00673         // annealing loop, stop when T<Tmin  (i.e. beta>1/Tmin)
00674         while (beta < betamax_) {
00675 
00676 
00677                 if(useTc_){
00678                   update(beta, tks,y, false, rho0);
00679                   while(merge(y,beta)){update(beta, tks,y, false, rho0);}
00680                   split(beta, tks,y);
00681                   beta=beta/coolingFactor_;
00682                 }else{
00683                   beta=beta/coolingFactor_;
00684                   splitAll(y);
00685                 }
00686 
00687                 // make sure we are not too far from equilibrium before cooling further
00688                 niter = 0;
00689                 while ((update(beta, tks, y, false, rho0) > 1.e-6) && (niter++
00690                                 < maxIterations_)) {
00691                 }
00692 
00693         }
00694 
00695 
00696         if(useTc_){
00697           // last round of splitting, make sure no critical clusters are left
00698           update(beta, tks,y, false, rho0);// make sure Tc is up-to-date
00699           while(merge(y,beta)){update(beta, tks,y, false, rho0);}
00700           unsigned int ntry=0;
00701           while( split(beta, tks,y) && (ntry++<10) ){
00702             niter=0; 
00703             while((update(beta, tks,y, false, rho0)>1.e-6)  && (niter++ < maxIterations_)){}
00704             merge(y,beta);
00705             update(beta, tks,y, false, rho0);
00706           }
00707         }else{
00708           // merge collapsed clusters 
00709           while(merge(y,beta)){update(beta, tks,y, false, rho0);}  
00710           //while(merge(y)){}   original code
00711         }
00712  
00713         if (verbose_) {
00714                 LogDebug("DAClusterizerinZ_vectorized")  << "dump after 1st merging " << endl;
00715                 dump(beta, y, tks, 2);
00716         }
00717 
00718         // switch on outlier rejection
00719         rho0 = 1. / nt;
00720         
00721         // auto-vectorized
00722         for (unsigned int k = 0; k < y.GetSize(); k++) {
00723                 y._pk[k] = 1.;
00724         } // democratic
00725         niter = 0;
00726         while ((update(beta, tks, y, true, rho0) > 1.e-8) && (niter++
00727                         < maxIterations_)) {
00728         }
00729         if (verbose_) {
00730                 LogDebug("DAClusterizerinZ_vectorized")  << "rho0=" << rho0 << " niter=" << niter << endl;
00731                 dump(beta, y, tks, 2);
00732         }
00733 
00734         // merge again  (some cluster split by outliers collapse here)
00735         while (merge(y)) {
00736         }
00737         if (verbose_) {
00738                 LogDebug("DAClusterizerinZ_vectorized")  << "dump after 2nd merging " << endl;
00739                 dump(beta, y, tks, 2);
00740         }
00741 
00742         // continue from freeze-out to Tstop (=1) without splitting, eliminate insignificant vertices
00743         while (beta <= betastop_) {
00744                 while (purge(y, tks, rho0, beta)) {
00745                         niter = 0;
00746                         while ((update(beta, tks, y, true, rho0) > 1.e-6) && (niter++
00747                                         < maxIterations_)) {
00748                         }
00749                 }
00750                 beta /= coolingFactor_;
00751                 niter = 0;
00752                 while ((update(beta, tks, y, true, rho0) > 1.e-6) && (niter++
00753                                 < maxIterations_)) {
00754                 }
00755         }
00756 
00757         if (verbose_) {
00758                 LogDebug("DAClusterizerinZ_vectorized")  << "Final result, rho0=" << rho0 << endl;
00759                 dump(beta, y, tks, 2);
00760         }
00761 
00762         // select significant tracks and use a TransientVertex as a container
00763         GlobalError dummyError;
00764 
00765         // ensure correct normalization of probabilities, should makes double assginment reasonably impossible
00766         const unsigned int nv = y.GetSize();
00767 
00768         for (unsigned int i = 0; i < nt; i++) {
00769                 tks._Z_sum[i] = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_);
00770 
00771                 for (unsigned int k = 0; k < nv; k++) {
00772                         tks._Z_sum[i] += y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k],
00773                                         tks._dz2[i]));
00774                 }
00775         }
00776 
00777         for (unsigned int k = 0; k < nv; k++) {
00778 
00779                 GlobalPoint pos(0, 0, y._z[k]);
00780 
00781                 vector<reco::TransientTrack> vertexTracks;
00782                 for (unsigned int i = 0; i < nt; i++) {
00783                         if (tks._Z_sum[i] > 0) {
00784 
00785                                 double p = y._pk[k] * local_exp(-beta * Eik(tks._z[i], y._z[k],
00786                                                 tks._dz2[i])) / tks._Z_sum[i];
00787                                 if ((tks._pi[i] > 0) && (p > 0.5)) {
00788 
00789                                         vertexTracks.push_back(*(tks.tt[i]));
00790                                         tks._Z_sum[i] = 0;
00791                                 } // setting Z=0 excludes double assignment
00792                         }
00793                 }
00794                 TransientVertex v(pos, dummyError, vertexTracks, 0);
00795                 clusters.push_back(v);
00796         }
00797 
00798         return clusters;
00799 
00800 }
00801 
00802 vector<vector<reco::TransientTrack> > DAClusterizerInZ_vect::clusterize(
00803                 const vector<reco::TransientTrack> & tracks) const {
00804                   
00805         if (verbose_) {
00806                 std::cout  << "###################################################" << endl;
00807                 std::cout  << "# vectorized DAClusterizerInZ_vect::clusterize   nt=" << tracks.size() << endl;
00808                 std::cout  << "###################################################" << endl;
00809         }
00810 
00811         vector<vector<reco::TransientTrack> > clusters;
00812         vector<TransientVertex> pv = vertices(tracks);
00813 
00814         if (verbose_) {
00815                 LogDebug("DAClusterizerinZ_vectorized")  << "# DAClusterizerInZ::clusterize   pv.size=" << pv.size()
00816                                 << endl;
00817         }
00818         if (pv.size() == 0) {
00819                 return clusters;
00820         }
00821 
00822         // fill into clusters and merge
00823         vector<reco::TransientTrack> aCluster = pv.begin()->originalTracks();
00824 
00825         for (vector<TransientVertex>::iterator k = pv.begin() + 1; k != pv.end(); k++) {
00826                 if ( fabs(k->position().z() - (k - 1)->position().z()) > (2 * vertexSize_)) {
00827                         // close a cluster
00828                         clusters.push_back(aCluster);
00829                         aCluster.clear();
00830                 }
00831                 for (unsigned int i = 0; i < k->originalTracks().size(); i++) {
00832                         aCluster.push_back(k->originalTracks().at(i));
00833                 }
00834 
00835         }
00836         clusters.push_back(aCluster);
00837 
00838         return clusters;
00839 
00840 }
00841