CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZ.cc

Go to the documentation of this file.
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   // prepare track data for clustering
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     //  get the beam-spot
00030     reco::BeamSpot beamspot=(it->stateAtBeamLine()).beamSpot();
00031     t.dz2= pow((*it).track().dzError(),2)          // track errror
00032       + (pow(beamspot.BeamWidthX(),2)+pow(beamspot.BeamWidthY(),2))/pow(tantheta,2)  // beam-width induced
00033       + pow(vertexSize_,2);                        // intrinsic vertex size, safer for outliers and short lived decays
00034     Measurement1D IP=(*it).stateAtBeamLine().transverseImpactParameter();// error constains beamspot
00035     t.pi=1./(1.+exp(pow(IP.value()/IP.error(),2)-pow(3.,2)));  // reduce weight for high ip tracks  
00036     t.tt=&(*it);
00037     t.Z=1.;
00038     tks.push_back(t);
00039   }
00040   return tks;
00041 }
00042 
00043 
00044 
00045 
00046 // double DAClusterizerInZ::pik(const double beta, const track_t & t, const vertex_t &k )const{
00047 //   //note: t.Z = sum_k exp(-beta*Eik ) is assumed to be valid, this is done in updateAndFit
00048 //   //      then we have sum_k pik(beta, t, k) = 1
00049 //   //      the last call of updateAndFit must have been made with the same temperature ! 
00050 //   //      at low T, the vertex position must be quite accurate (iterated)
00051 //   double Eik=pow(t.z-k.z,2)/t.dz2;
00052 //   double o= k.pk*exp(-beta*Eik); 
00053 //   if(isnan(o)){ cout << "NaN in pik " << k.pk << " " << beta << " " << Eik << " Z=" << t.Z << "   beta Eik=" << beta*Eik << endl;}
00054 //   if (t.Z>0) return k.pk*exp(-beta*Eik)/t.Z; 
00055 //   return 0;
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   //update weights and vertex positions
00072   // mass constrained annealing without noise
00073   // returns the squared sum of changes of vertex positions
00074 
00075   //cout << endl << endl << "DAClusterizerInZ::update" << endl;
00076   unsigned int nt=tks.size();
00077 
00078   //initialize sums
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   // loop over tracks
00086   for(unsigned int i=0; i<nt; i++){
00087 
00088     // update pik and Zi
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));// cache exponential for one track at a time
00092       Zi += k->pk * k->ei;
00093     }
00094     tks[i].Z=Zi;
00095 
00096     // normalization for pk
00097     if (tks[i].Z>0){
00098       sumpi += tks[i].pi;
00099       // accumulate weighted z and weights for vertex update
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   } // end of track loop
00112 
00113 
00114   // now update z and pk
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   // return how much the prototypes moved
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   // MVF style, no more vertex weights, update tracks weights and vertex positions, with noise 
00144   // returns the squared sum of changes of vertex positions
00145 
00146   unsigned int nt=tks.size();
00147 
00148   //initialize sums
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   // loop over tracks
00155   for(unsigned int i=0; i<nt; i++){
00156 
00157     // update pik and Zi
00158     double Zi=rho0*exp(-beta*mu0_*mu0_);// cut-off
00159     for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
00160       k->ei=exp(-beta*Eik(tks[i],*k));// cache exponential for one track at a time
00161       Zi += k->pk * k->ei;
00162     }
00163     tks[i].Z=Zi;
00164 
00165     // normalization
00166     if (tks[i].Z>0){
00167       //double pi0=exp(-beta*mu0_*mu0_)/tks[i].Z;
00168       //sumpi += tks[i].pi*(1.-pi0*rho0);  // exclude rho0 from the normalization of the cluster "mass"
00169        // accumulate weighted z and weights for vertex update
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     //else{  sumpi += tks[i].pi;  }
00178 
00179 
00180   } // end of track loop
00181 
00182 
00183   // now update z and pk
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     //k->pk *= k->se / sumpi; 
00196   }
00197 
00198   // return how much the prototypes moved
00199   return delta;
00200 }
00201 
00202 
00203 
00204 
00205 
00206 bool DAClusterizerInZ::merge(vector<vertex_t> & y, int nt)const{
00207   // merge clusters that collapsed or never separated, return true if vertices were merged, false otherwise
00208 
00209  if(y.size()<2)  return false;
00210 
00211   for(vector<vertex_t>::iterator k=y.begin(); (k+1)!=y.end(); k++){
00212     //if ((k+1)->z - k->z<1.e-2){  // note, no fabs here, maintains z-ordering  (with split()+merge() at every temperature)
00213     if (fabs((k+1)->z - k->z)<1.e-2){  // with fabs if only called after freeze-out (splitAll() at highter T)
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   // eliminate clusters with only one significant/unique track
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         //double p=pik(beta,tks[i],*k);
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     //rho0+=k0->pk;
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;  // max Tc for beta=0
00272   // estimate critical temperature from beta=0 (T=inf)
00273   unsigned int nt=tks.size();
00274 
00275   for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
00276 
00277     // vertex fit at T=inf 
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     // estimate Tcrit, eventually do this in the same loop
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;  // the critical temperature of this vertex
00296     if(Tc>T0) T0=Tc;
00297   }// vertex loop (normally there should be only one vertex at beta=0)
00298   
00299   if (T0>1./betamax){
00300     return betamax/pow(coolingFactor_, int(log(T0*betamax)/log(coolingFactor_))-1 );
00301   }else{
00302     // ensure at least one annealing step
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;      // split all single vertices by 10 um 
00317   double zsep=2*epsilon;    // split vertices that are isolated by at least zsep (vertices that haven't collapsed)
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       // isolated prototype, split
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   }// vertex loop
00339   
00340   y=y1;
00341 }
00342  
00343 
00344 DAClusterizerInZ::DAClusterizerInZ(const edm::ParameterSet& conf) 
00345 {
00346   // some defaults to avoid uninitialized variables
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;  // 0.5 mm
00353   mu0_=4.0;   // Adaptive Fitter uses 3.0 but that appears to be a bit tight here sometimes
00354 
00355   // configure
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   // copy and sort for nicer printout
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(); // see DataFormats/TrackReco/interface/HitPattern.h
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           //double p=pik(beta,tks[i],*k);
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;  // start with no outlier rejection
00450 
00451   vector< TransientVertex > clusters;
00452   if (tks.empty()) return clusters;
00453 
00454   vector<vertex_t> y; // the vertex prototypes
00455 
00456   // initialize:single vertex at infinite temperature
00457   vertex_t vstart;
00458   vstart.z=0.;
00459   vstart.pk=1.;
00460   y.push_back(vstart);
00461   int niter=0;      // number of iterations
00462   
00463 
00464   // estimate first critical temperature
00465   double beta=beta0(betamax_, tks, y);
00466   niter=0; while((update(beta, tks,y)>1.e-6)  && (niter++ < maxIterations_)){ }
00467 
00468   
00469  // annealing loop, stop when T<Tmin  (i.e. beta>1/Tmin)
00470   while(beta<betamax_){ 
00471 
00472     beta=beta/coolingFactor_;
00473     splitAll(tks,y);
00474 
00475     // make sure we are not too far from equilibrium before cooling further
00476     niter=0; while((update(beta, tks,y)>1.e-6)  && (niter++ < maxIterations_)){ }
00477 
00478   }
00479 
00480 
00481   // merge collapsed clusters 
00482   while(merge(y,tks.size())){} 
00483   if(verbose_  ){ cout << "dump after 1st merging " << endl;  dump(beta,y,tks,2);}
00484 
00485 
00486 
00487   // switch on outlier rejection
00488   //rho0=exp(beta*mu0_*mu0_)/nt; if(rho0>0.1) rho0=0.1;
00489   rho0=1./nt; for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){ k->pk =1.; }  // democratic
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   // merge again  (some cluster split by outliers collapse here)
00495   while(merge(y,tks.size())){}  
00496   if(verbose_  ){ cout << "dump after 2nd merging " << endl;  dump(beta,y,tks,2);}
00497 
00498 
00499   // continue from freeze-out to Tstop (=1) without splitting, eliminate insignificant vertices
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   // select significant tracks and use a TransientVertex as a container
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   // fill into clusters and merge
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       // close a cluster
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