00001 #include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ.h"
00002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00003 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
00004 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
00005
00006
00007 using namespace std;
00008
00009
00010 namespace {
00011
00012 bool recTrackLessZ1(const DAClusterizerInZ::track_t & tk1,
00013 const DAClusterizerInZ::track_t & tk2)
00014 {
00015 return tk1.z < tk2.z;
00016 }
00017 }
00018
00019
00020 vector<DAClusterizerInZ::track_t> DAClusterizerInZ::fill(
00021 const vector<reco::TransientTrack> & tracks
00022 )const{
00023
00024 vector<track_t> tks;
00025 for(vector<reco::TransientTrack>::const_iterator it=tracks.begin(); it!=tracks.end(); it++){
00026 track_t t;
00027 t.z=((*it).stateAtBeamLine().trackStateAtPCA()).position().z();
00028 double tantheta=tan(((*it).stateAtBeamLine().trackStateAtPCA()).momentum().theta());
00029
00030 reco::BeamSpot beamspot=(it->stateAtBeamLine()).beamSpot();
00031 t.dz2= pow((*it).track().dzError(),2)
00032 + (pow(beamspot.BeamWidthX(),2)+pow(beamspot.BeamWidthY(),2))/pow(tantheta,2)
00033 + pow(vertexSize_,2);
00034 Measurement1D IP=(*it).stateAtBeamLine().transverseImpactParameter();
00035 t.pi=1./(1.+exp(pow(IP.value()/IP.error(),2)-pow(3.,2)));
00036 t.tt=&(*it);
00037 t.Z=1.;
00038 tks.push_back(t);
00039 }
00040 return tks;
00041 }
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 double DAClusterizerInZ::Eik(const track_t & t, const vertex_t &k )const{
00060 return pow(t.z-k.z,2)/t.dz2;
00061 }
00062
00063
00064
00065
00066 double DAClusterizerInZ::update(
00067 double beta,
00068 vector<track_t> & tks,
00069 vector<vertex_t> & y
00070 )const{
00071
00072
00073
00074
00075
00076 unsigned int nt=tks.size();
00077
00078
00079 double sumpi=0;
00080 for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
00081 k->se=0; k->sw=0; k->swz=0;
00082 }
00083
00084
00085
00086 for(unsigned int i=0; i<nt; i++){
00087
00088
00089 double Zi=0.;
00090 for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
00091 k->ei=exp(-beta*Eik(tks[i],*k));
00092 Zi += k->pk * k->ei;
00093 }
00094 tks[i].Z=Zi;
00095
00096
00097 if (tks[i].Z>0){
00098 sumpi += tks[i].pi;
00099
00100 for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
00101 k->se += tks[i].pi* k->ei / Zi;
00102 double w = k->pk * tks[i].pi* k->ei / Zi / tks[i].dz2;
00103 k->sw += w;
00104 k->swz += w * tks[i].z;
00105 }
00106 }else{
00107 sumpi += tks[i].pi;
00108 }
00109
00110
00111 }
00112
00113
00114
00115 double delta=0;
00116 for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
00117 if ( k->sw > 0){
00118 double znew=k->swz/k->sw;
00119 delta+=pow(k->z-znew,2);
00120 k->z=znew;
00121 }else{
00122 edm::LogInfo("sumw") << "invalid sum of weights in fit: " << k->sw << endl;
00123 if(verbose_){cout << " a cluster melted away ? pk=" << k->pk << " sumw=" << k->sw << endl;}
00124 }
00125
00126 k->pk = k->pk * k->se / sumpi;
00127 }
00128
00129
00130 return delta;
00131 }
00132
00133
00134
00135
00136
00137 double DAClusterizerInZ::update(
00138 double beta,
00139 vector<track_t> & tks,
00140 vector<vertex_t> & y,
00141 double & rho0
00142 )const{
00143
00144
00145
00146 unsigned int nt=tks.size();
00147
00148
00149 for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
00150 k->se=0; k->sw=0; k->swz=0;
00151 }
00152
00153
00154
00155 for(unsigned int i=0; i<nt; i++){
00156
00157
00158 double Zi=rho0*exp(-beta*mu0_*mu0_);
00159 for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
00160 k->ei=exp(-beta*Eik(tks[i],*k));
00161 Zi += k->pk * k->ei;
00162 }
00163 tks[i].Z=Zi;
00164
00165
00166 if (tks[i].Z>0){
00167
00168
00169
00170 for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
00171 k->se += tks[i].pi* k->ei / Zi;
00172 double w = k->pk * tks[i].pi* k->ei / Zi / tks[i].dz2;
00173 k->sw += w;
00174 k->swz += w * tks[i].z;
00175 }
00176 }
00177
00178
00179
00180 }
00181
00182
00183
00184 double delta=0;
00185 for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
00186 if ( k->sw > 0){
00187 double znew=k->swz/k->sw;
00188 delta+=pow(k->z-znew,2);
00189 k->z=znew;
00190 }else{
00191 edm::LogInfo("sumw") << "invalid sum of weights in fit: " << k->sw << endl;
00192 if(verbose_){cout << " a cluster melted away ? pk=" << k->pk << " sumw=" << k->sw << endl;}
00193 }
00194
00195
00196 }
00197
00198
00199 return delta;
00200 }
00201
00202
00203
00204
00205
00206 bool DAClusterizerInZ::merge(vector<vertex_t> & y, int nt)const{
00207
00208
00209 if(y.size()<2) return false;
00210
00211 for(vector<vertex_t>::iterator k=y.begin(); (k+1)!=y.end(); k++){
00212
00213 if (fabs((k+1)->z - k->z)<1.e-2){
00214 k->pk += (k+1)->pk;
00215 k->z=0.5*(k->z+(k+1)->z);
00216 y.erase(k+1);
00217 return true;
00218 }
00219 }
00220
00221 return false;
00222 }
00223
00224
00225
00226 bool DAClusterizerInZ::purge(vector<vertex_t> & y, vector<track_t> & tks, double & rho0, const double beta)const{
00227
00228 if(y.size()<2) return false;
00229
00230 unsigned int nt=tks.size();
00231 double sumpmin=nt;
00232 vector<vertex_t>::iterator k0=y.end();
00233 for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
00234 int nUnique=0;
00235 double sump=0;
00236 double pmax=k->pk/(k->pk+rho0*exp(-beta*mu0_*mu0_));
00237 for(unsigned int i=0; i<nt; i++){
00238 if(tks[i].Z>0){
00239
00240 double p=k->pk * exp(-beta*Eik(tks[i],*k)) / tks[i].Z;
00241 sump+=p;
00242 if( (p > 0.9*pmax) && (tks[i].pi>0) ){ nUnique++; }
00243 }
00244 }
00245
00246 if((nUnique<2)&&(sump<sumpmin)){
00247 sumpmin=sump;
00248 k0=k;
00249 }
00250 }
00251
00252 if(k0!=y.end()){
00253 if(verbose_){cout << "eliminating prototype at " << k0->z << " with sump=" << sumpmin << endl;}
00254
00255 y.erase(k0);
00256 return true;
00257 }else{
00258 return false;
00259 }
00260 }
00261
00262
00263
00264
00265 double DAClusterizerInZ::beta0(
00266 double betamax,
00267 vector<track_t> & tks,
00268 vector<vertex_t> & y
00269 )const{
00270
00271 double T0=0;
00272
00273 unsigned int nt=tks.size();
00274
00275 for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
00276
00277
00278 double sumwz=0;
00279 double sumw=0;
00280 for(unsigned int i=0; i<nt; i++){
00281 double w=tks[i].pi/tks[i].dz2;
00282 sumwz+=w*tks[i].z;
00283 sumw +=w;
00284 }
00285 k->z=sumwz/sumw;
00286
00287
00288 double a=0, b=0;
00289 for(unsigned int i=0; i<nt; i++){
00290 double dx=tks[i].z-(k->z);
00291 double w=tks[i].pi/tks[i].dz2;
00292 a+=w*pow(dx,2)/tks[i].dz2;
00293 b+=w;
00294 }
00295 double Tc= 2.*a/b;
00296 if(Tc>T0) T0=Tc;
00297 }
00298
00299 if (T0>1./betamax){
00300 return betamax/pow(coolingFactor_, int(log(T0*betamax)/log(coolingFactor_))-1 );
00301 }else{
00302
00303 return betamax/coolingFactor_;
00304 }
00305 }
00306
00307
00308
00309
00310 void DAClusterizerInZ::splitAll(
00311 vector<track_t> & tks,
00312 vector<vertex_t> & y
00313 )const{
00314
00315
00316 double epsilon=1e-3;
00317 double zsep=2*epsilon;
00318 vector<vertex_t> y1;
00319
00320 for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
00321 if ( ( (k==y.begin())|| (k-1)->z < k->z - zsep) && (((k+1)==y.end() )|| (k+1)->z > k->z + zsep)) {
00322
00323 vertex_t vnew;
00324 vnew.z = k->z -epsilon;
00325 (*k).z = k->z+epsilon;
00326 vnew.pk= 0.5* (*k).pk;
00327 (*k).pk= 0.5* (*k).pk;
00328 y1.push_back(vnew);
00329 y1.push_back(*k);
00330
00331 }else if( y1.empty() || (y1.back().z < k->z -zsep)){
00332 y1.push_back(*k);
00333 }else{
00334 y1.back().z -=epsilon;
00335 k->z+=epsilon;
00336 y1.push_back(*k);
00337 }
00338 }
00339
00340 y=y1;
00341 }
00342
00343
00344 DAClusterizerInZ::DAClusterizerInZ(const edm::ParameterSet& conf)
00345 {
00346
00347 verbose_= conf.getUntrackedParameter<bool>("verbose", false);
00348 betamax_=0.1;
00349 betastop_ =1.0;
00350 coolingFactor_=0.8;
00351 maxIterations_=100;
00352 vertexSize_=0.05;
00353 mu0_=4.0;
00354
00355
00356
00357 double Tmin = conf.getParameter<double>("Tmin");
00358 vertexSize_ = conf.getParameter<double>("vertexSize");
00359 coolingFactor_ = conf.getParameter<double>("coolingFactor");
00360 maxIterations_=100;
00361 if (Tmin==0){
00362 cout << "DAClusterizerInZ: invalid Tmin" << Tmin << " reset do default " << 1./betamax_ << endl;
00363 }else{
00364 betamax_ = 1./Tmin;
00365 }
00366
00367 }
00368
00369
00370 void DAClusterizerInZ::dump(const double beta, const vector<vertex_t> & y, const vector<track_t> & tks0, int verbosity)const{
00371
00372
00373 vector<track_t> tks;
00374 for(vector<track_t>::const_iterator t=tks0.begin(); t!=tks0.end(); t++){tks.push_back(*t); }
00375 stable_sort(tks.begin(), tks.end(), recTrackLessZ1);
00376
00377 cout << "-----DAClusterizerInZ::dump ----" << endl;
00378 cout << "beta=" << beta << " betamax= " << betamax_ << endl;
00379 cout << " z= ";
00380 cout.precision(4);
00381 for(vector<vertex_t>::const_iterator k=y.begin(); k!=y.end(); k++){
00382 cout << setw(8) << fixed << k->z ;
00383 }
00384 cout << endl << "T=" << setw(15) << 1./beta <<" Tc= ";
00385 cout << endl << " pk=";
00386 double sumpk=0;
00387 for(vector<vertex_t>::const_iterator k=y.begin(); k!=y.end(); k++){
00388 cout << setw(8) << setprecision(3) << fixed << k->pk;
00389 sumpk+=k->pk;
00390 }
00391 cout << endl;
00392
00393 if(verbosity>0){
00394 double E=0, F=0;
00395 cout << endl;
00396 cout << "---- z +/- dz ip +/-dip pt phi eta weights ----" << endl;
00397 cout.precision(4);
00398 for(unsigned int i=0; i<tks.size(); i++){
00399 if (tks[i].Z>0){ F-=log(tks[i].Z)/beta;}
00400 double tz= tks[i].z;
00401 cout << setw (3)<< i << ")" << setw (8) << fixed << setprecision(4)<< tz << " +/-" << setw (6)<< sqrt(tks[i].dz2);
00402
00403 if(tks[i].tt->track().quality(reco::TrackBase::highPurity)){ cout << " *";}else{cout <<" ";}
00404 if(tks[i].tt->track().hitPattern().hasValidHitInFirstPixelBarrel()){cout <<"+";}else{cout << "-";}
00405 cout << setw(1) << tks[i].tt->track().hitPattern().pixelBarrelLayersWithMeasurement();
00406 cout << setw(1) << tks[i].tt->track().hitPattern().pixelEndcapLayersWithMeasurement();
00407 cout << setw(1) << hex << tks[i].tt->track().hitPattern().trackerLayersWithMeasurement()-tks[i].tt->track().hitPattern().pixelLayersWithMeasurement() <<dec;
00408 cout << "=" << setw(1)<<hex <<tks[i].tt->track().trackerExpectedHitsOuter().numberOfHits() << dec;
00409
00410 Measurement1D IP=tks[i].tt->stateAtBeamLine().transverseImpactParameter();
00411 cout << setw (8) << IP.value() << "+/-" << setw (6) << IP.error();
00412 cout << " " << setw(6) << setprecision(2) << tks[i].tt->track().pt()*tks[i].tt->track().charge();
00413 cout << " " << setw(5) << setprecision(2) << tks[i].tt->track().phi()
00414 << " " << setw(5) << setprecision(2) << tks[i].tt->track().eta() ;
00415
00416 double sump=0.;
00417 for(vector<vertex_t>::const_iterator k=y.begin(); k!=y.end(); k++){
00418 if((tks[i].pi>0)&&(tks[i].Z>0)){
00419
00420 double p=k->pk * exp(-beta*Eik(tks[i],*k)) / tks[i].Z;
00421 if( p > 0.0001){
00422 cout << setw (8) << setprecision(3) << p;
00423 }else{
00424 cout << " . ";
00425 }
00426 E+=p*Eik(tks[i],*k);
00427 sump+=p;
00428 }else{
00429 cout << " ";
00430 }
00431 }
00432 cout << endl;
00433 }
00434 cout << endl << "T=" << 1/beta << " E=" << E << " n="<< y.size() << " F= " << F << endl << "----------" << endl;
00435 }
00436 }
00437
00438
00439
00440
00441
00442 vector< TransientVertex >
00443 DAClusterizerInZ::vertices(const vector<reco::TransientTrack> & tracks, const int verbosity)
00444 const
00445 {
00446
00447 vector<track_t> tks=fill(tracks);
00448 unsigned int nt=tracks.size();
00449 double rho0=0.0;
00450
00451 vector< TransientVertex > clusters;
00452 if (tks.empty()) return clusters;
00453
00454 vector<vertex_t> y;
00455
00456
00457 vertex_t vstart;
00458 vstart.z=0.;
00459 vstart.pk=1.;
00460 y.push_back(vstart);
00461 int niter=0;
00462
00463
00464
00465 double beta=beta0(betamax_, tks, y);
00466 niter=0; while((update(beta, tks,y)>1.e-6) && (niter++ < maxIterations_)){ }
00467
00468
00469
00470 while(beta<betamax_){
00471
00472 beta=beta/coolingFactor_;
00473 splitAll(tks,y);
00474
00475
00476 niter=0; while((update(beta, tks,y)>1.e-6) && (niter++ < maxIterations_)){ }
00477
00478 }
00479
00480
00481
00482 while(merge(y,tks.size())){}
00483 if(verbose_ ){ cout << "dump after 1st merging " << endl; dump(beta,y,tks,2);}
00484
00485
00486
00487
00488
00489 rho0=1./nt; for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){ k->pk =1.; }
00490 niter=0; while((update(beta, tks,y,rho0) > 1.e-8) && (niter++ < maxIterations_)){ }
00491 if(verbose_ ){ cout << "rho0=" << rho0 << " niter=" << niter << endl; dump(beta,y,tks,2);}
00492
00493
00494
00495 while(merge(y,tks.size())){}
00496 if(verbose_ ){ cout << "dump after 2nd merging " << endl; dump(beta,y,tks,2);}
00497
00498
00499
00500 while(beta<=betastop_){
00501 while(purge(y,tks,rho0, beta)){
00502 niter=0; while((update(beta, tks,y,rho0) > 1.e-6) && (niter++ < maxIterations_)){ }
00503 }
00504 beta/=coolingFactor_;
00505 niter=0; while((update(beta, tks,y,rho0) > 1.e-6) && (niter++ < maxIterations_)){ }
00506 }
00507
00508 if(verbose_){
00509 cout << "Final result, rho0=" << rho0 << endl;
00510 dump(beta,y,tks,2);
00511 }
00512
00513
00514
00515 GlobalError dummyError;
00516 for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
00517 GlobalPoint pos(0, 0, k->z);
00518 vector< reco::TransientTrack > vertexTracks;
00519 for(unsigned int i=0; i<nt; i++){
00520 if(tks[i].Z>0){
00521 double p=k->pk * exp(-beta*Eik(tks[i],*k)) / tks[i].Z;
00522 if( (tks[i].pi>0) && ( p > 0.5 ) ){ vertexTracks.push_back(*(tks[i].tt));}
00523 }
00524 }
00525 TransientVertex v(pos, dummyError, vertexTracks, 0);
00526 clusters.push_back(v);
00527 }
00528
00529 return clusters;
00530
00531 }
00532
00533
00534
00535
00536
00537 vector< vector<reco::TransientTrack> >
00538 DAClusterizerInZ::clusterize(const vector<reco::TransientTrack> & tracks)
00539 const
00540 {
00541 if(verbose_) {
00542 cout << "###################################################" << endl;
00543 cout << "# DAClusterizerInZ::clusterize nt="<<tracks.size() << endl;
00544 cout << "###################################################" << endl;
00545 }
00546
00547 vector< vector<reco::TransientTrack> > clusters;
00548 vector< TransientVertex > pv=vertices(tracks);
00549
00550 if(verbose_){ cout << "# DAClusterizerInZ::clusterize pv.size="<<pv.size() << endl; }
00551 if (pv.size()==0){ return clusters;}
00552
00553
00554
00555 vector< reco::TransientTrack> aCluster=pv.begin()->originalTracks();
00556
00557 for(vector<TransientVertex>::iterator k=pv.begin()+1; k!=pv.end(); k++){
00558 if ( k->position().z() - (k-1)->position().z()> (2*vertexSize_) ){
00559
00560 clusters.push_back(aCluster);
00561 aCluster.clear();
00562 }
00563 for(unsigned int i=0; i<k->originalTracks().size(); i++){ aCluster.push_back( k->originalTracks().at(i)); }
00564
00565 }
00566 clusters.push_back(aCluster);
00567
00568
00569 return clusters;
00570
00571 }
00572