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
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;
00026 dzCutOff_ = 4.0;
00027
00028
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
00062 for(int i=0; i!=arg_arr_size; ++i ) arg_out[i]=vdt::fast_exp(arg_inp[i]);
00063
00064 }
00065 else
00066 for(int i=0; i!=arg_arr_size; ++i ) arg_out[i]=std::exp(arg_inp[i]);
00067 }
00068
00069
00070 DAClusterizerInZ_vect::track_t DAClusterizerInZ_vect::fill(const vector<
00071 reco::TransientTrack> & tracks) const {
00072
00073
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
00084 reco::BeamSpot beamspot = (it->stateAtBeamLine()).beamSpot();
00085 double t_dz2 =
00086 pow((*it).track().dzError(), 2)
00087 + (pow(beamspot.BeamWidthX()*cos(phi),2)+pow(beamspot.BeamWidthY()*sin(phi),2))/pow(tantheta,2)
00088 + pow(vertexSize_, 2);
00089 if (d0CutOff_ > 0) {
00090 Measurement1D IP =
00091 (*it).stateAtBeamLine().transverseImpactParameter();
00092 t_pi = 1. / (1. + local_exp(pow(IP.value() / IP.error(), 2) - pow(
00093 d0CutOff_, 2)));
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
00120
00121
00122
00123 const unsigned int nt = gtracks.GetSize();
00124 const unsigned int nv = gvertices.GetSize();
00125
00126
00127 double sumpi = 0;
00128
00129
00130 double delta = 0;
00131
00132 unsigned int itrack;
00133 unsigned int ivertex;
00134
00135
00136 double Z_init = 0;
00137
00138
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
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
00157
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
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
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
00200 if ( useRho0 )
00201 {
00202 Z_init = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_);
00203 }
00204
00205
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
00210
00211 gtracks._Z_sum[itrack] = kernel_add_Z(gvertices);
00212
00213
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
00223 auto kernel_calc_z = [ &delta, &sumpi, nv, this, useRho0 ] (vertex_t & vertices )
00224 {
00225
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
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
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
00256 return delta;
00257 }
00258
00259
00260
00261 bool DAClusterizerInZ_vect::merge(vertex_t & y, double & beta)const{
00262
00263
00264
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
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
00309 if (fabs(y._z[k + 1] - y._z[k]) < 1.e-2) {
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
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
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
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;
00381
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
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
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;
00406 if (Tc > T0)
00407 T0 = Tc;
00408 }
00409
00410 if (T0 > 1. / betamax) {
00411 return betamax / pow(coolingFactor_, int(log(T0 * betamax) / log(
00412 coolingFactor_)) - 1);
00413 } else {
00414
00415 return betamax / coolingFactor_;
00416 }
00417 }
00418
00419
00420 bool DAClusterizerInZ_vect::split(const double beta, track_t &tks, vertex_t & y ) const{
00421
00422
00423
00424
00425 double epsilon=1e-3;
00426 unsigned int nv = y.GetSize();
00427
00428
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
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
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
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
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;
00495 double zsep = 2 * epsilon;
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
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 }
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();
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
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;
00648
00649 vector<TransientVertex> clusters;
00650 if (tks.GetSize() == 0)
00651 return clusters;
00652
00653 vertex_t y;
00654
00655
00656 y.AddItem( 0, 1.0);
00657
00658 int niter = 0;
00659
00660
00661
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
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
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
00698 update(beta, tks,y, false, rho0);
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
00709 while(merge(y,beta)){update(beta, tks,y, false, rho0);}
00710
00711 }
00712
00713 if (verbose_) {
00714 LogDebug("DAClusterizerinZ_vectorized") << "dump after 1st merging " << endl;
00715 dump(beta, y, tks, 2);
00716 }
00717
00718
00719 rho0 = 1. / nt;
00720
00721
00722 for (unsigned int k = 0; k < y.GetSize(); k++) {
00723 y._pk[k] = 1.;
00724 }
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
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
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
00763 GlobalError dummyError;
00764
00765
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 }
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
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
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