Fork me on GitHub

Changeset 4154bbd in git


Ignore:
Timestamp:
Sep 1, 2016, 12:56:40 PM (8 years ago)
Author:
Pavel Demin <pavel.demin@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
61569e0
Parents:
95b4e9f
Message:

adapt vertexing modules to ROOT 5.34 and GCC 4.4

Files:
4 edited

Legend:

Unmodified
Added
Removed

  • Makefile

    r95b4e9f r4154bbd  
    16691669
    16701670modules/VertexFinderDA4D.h: \
    1671         classes/DelphesModule.h \
    1672         classes/DelphesClasses.h
     1671        classes/DelphesModule.h
    16731672        @touch $@
    16741673
    … …  
    19421941
    19431942modules/VertexSorter.h: \
    1944         classes/DelphesModule.h \
    1945         classes/DelphesClasses.h
     1943        classes/DelphesModule.h
    19461944        @touch $@
    19471945
    … …  
    19511949
    19521950modules/VertexFinder.h: \
    1953         classes/DelphesModule.h \
    1954         classes/DelphesClasses.h
     1951        classes/DelphesModule.h
    19551952        @touch $@
    19561953

  • modules/VertexFinderDA4D.cc

    r95b4e9f r4154bbd  
    4242static const Double_t c_light   = 2.99792458e+8 * m/s;
    4343
    44 
    45 namespace {
    46   bool recTrackLessZ1(const VertexFinderDA4D::track_t & tk1,
    47                       const VertexFinderDA4D::track_t & tk2)
    48   {
    49     return tk1.z < tk2.z;
    50   }
     44struct track_t
     45{
     46  double z;      // z-coordinate at point of closest approach to the beamline
     47  double t;      // t-coordinate at point of closest approach to the beamline
     48  double dz2;    // square of the error of z(pca)
     49  double dtz;    // covariance of z-t
     50  double dt2;    // square of the error of t(pca)
     51  Candidate *tt; // a pointer to the Candidate Track
     52  double Z;      // Z[i]   for DA clustering
     53  double pi;     // track weight
     54  double pt;
     55  double eta;
     56  double phi;
     57};
     58
     59struct vertex_t
     60{
     61  double z;
     62  double t;
     63  double pk;     // vertex weight for "constrained" clustering
     64  // --- temporary numbers, used during update
     65  double ei;
     66  double sw;
     67  double swz;
     68  double swt;
     69  double se;
     70  // ---for Tc
     71  double swE;
     72  double Tc;
     73};
     74
     75static bool split(double beta, std::vector<track_t> &tks, std::vector<vertex_t> &y);
     76static double update1(double beta, std::vector<track_t> &tks, std::vector<vertex_t> &y);
     77static double update2(double beta, std::vector<track_t> &tks, std::vector<vertex_t> &y, double &rho0, const double dzCutOff);
     78static void dump(const double beta, const std::vector<vertex_t> & y, const std::vector<track_t> & tks);
     79static bool merge(std::vector<vertex_t> &);
     80static bool merge(std::vector<vertex_t> &, double &);
     81static bool purge(std::vector<vertex_t> &, std::vector<track_t> & , double &, const double, const double);
     82static void splitAll(std::vector<vertex_t> &y);
     83static double beta0(const double betamax, std::vector<track_t> &tks, std::vector<vertex_t> &y, const double coolingFactor);
     84static double Eik(const track_t &t, const vertex_t &k);
     85
     86static bool recTrackLessZ1(const track_t & tk1, const track_t & tk2)
     87{
     88  return tk1.z < tk2.z;
    5189}
    5290
    … …  
    214252
    215253     if (fVerbose){
    216            std::cout << "x,y,z";
     254     std::cout << "x,y,z";
    217255       std::cout << ",t";
    218256       std::cout << "=" << candidate->Position.X()/10.0 <<" " << candidate->Position.Y()/10.0 << " " <<  candidate->Position.Z()/10.0;
    … …  
    242280//------------------------------------------------------------------------------
    243281
    244 
    245 void VertexFinderDA4D::clusterize(const TObjArray & tracks, TObjArray & clusters)
     282void VertexFinderDA4D::clusterize(const TObjArray &tracks, TObjArray &clusters)
    246283{
    247284  if(fVerbose) {
    … …  
    294331}
    295332
    296 
    297333//------------------------------------------------------------------------------
    298334
    … …  
    303339  vector< Candidate* > clusters;
    304340
    305   vector<track_t> tks = fill();
     341  vector<track_t> tks;
     342  track_t tr;
     343  Double_t z, dz, t, l, dt, d0, d0error;
     344
     345  // loop over input tracks
     346  fItInputArray->Reset();
     347  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
     348  {
     349    //TBC everything in cm
     350    z = candidate->DZ/10;
     351    tr.z = z;
     352    dz = candidate->ErrorDZ/10;
     353    tr.dz2 = dz*dz          // track error
     354    //TBC: beamspot size induced error, take 0 for now.
     355    // + (std::pow(beamspot.BeamWidthX()*cos(phi),2.)+std::pow(beamspot.BeamWidthY()*sin(phi),2.))/std::pow(tantheta,2.) // beam-width induced
     356    + fVertexSpaceSize*fVertexSpaceSize; // intrinsic vertex size, safer for outliers and short lived decays
     357
     358    // TBC: the time is in ns for now TBC
     359    //t = candidate->Position.T()/c_light;
     360    t = candidate->InitialPosition.T()/c_light;
     361    l = candidate->L/c_light;
     362    double pt = candidate->Momentum.Pt();
     363    double eta = candidate->Momentum.Eta();
     364    double phi = candidate->Momentum.Phi();
     365
     366    tr.pt = pt;
     367    tr.eta = eta;
     368    tr.phi = phi;
     369    tr.t = t; //
     370    tr.dtz = 0.;
     371    dt = candidate->ErrorT/c_light;
     372    tr.dt2 = dt*dt + fVertexTimeSize*fVertexTimeSize;   // the ~injected~ timing error plus a small minimum vertex size in time
     373    if(fD0CutOff>0)
     374    {
     375
     376      d0 = TMath::Abs(candidate->D0)/10.0;
     377      d0error = candidate->ErrorD0/10.0;
     378
     379      tr.pi=1./(1.+exp((d0*d0)/(d0error*d0error) - fD0CutOff*fD0CutOff));  // reduce weight for high ip tracks
     380
     381    }
     382    else
     383    {
     384      tr.pi=1.;
     385    }
     386    tr.tt=&(*candidate);
     387    tr.Z=1.;
     388
     389    // TBC now putting track selection here (> fPTMin)
     390    if(tr.pi > 1e-3 && tr.pt > fMinPT)
     391    {
     392      tks.push_back(tr);
     393    }
     394  }
    306395
    307396  //print out input tracks
    … …  
    341430
    342431  // estimate first critical temperature
    343   double beta=beta0(fBetaMax, tks, y);
    344   niter=0; while((update(beta, tks,y)>1.e-6)  && (niter++ < fMaxIterations)){ }
     432  double beta=beta0(fBetaMax, tks, y, fCoolingFactor);
     433  niter=0; while((update1(beta, tks,y)>1.e-6)  && (niter++ < fMaxIterations)){ }
    345434
    346435  // annealing loop, stop when T<Tmin  (i.e. beta>1/Tmin)
    … …  
    348437
    349438    if(fUseTc){
    350       update(beta, tks,y);
    351       while(merge(y,beta)){update(beta, tks,y);}
     439      update1(beta, tks,y);
     440      while(merge(y,beta)){update1(beta, tks,y);}
    352441      split(beta, tks,y);
    353442      beta=beta/fCoolingFactor;
    … …  
    358447
    359448   // make sure we are not too far from equilibrium before cooling further
    360    niter=0; while((update(beta, tks,y)>1.e-6)  && (niter++ < fMaxIterations)){ }
     449   niter=0; while((update1(beta, tks,y)>1.e-6)  && (niter++ < fMaxIterations)){ }
    361450
    362451  }
    … …  
    364453  if(fUseTc){
    365454    // last round of splitting, make sure no critical clusters are left
    366     update(beta, tks,y);
    367     while(merge(y,beta)){update(beta, tks,y);}
     455    update1(beta, tks,y);
     456    while(merge(y,beta)){update1(beta, tks,y);}
    368457    unsigned int ntry=0;
    369458    while( split(beta, tks,y) && (ntry++<10) ){
    370459      niter=0;
    371       while((update(beta, tks,y)>1.e-6)  && (niter++ < fMaxIterations)){}
     460      while((update1(beta, tks,y)>1.e-6)  && (niter++ < fMaxIterations)){}
    372461      merge(y,beta);
    373       update(beta, tks,y);
     462      update1(beta, tks,y);
    374463    }
    375464  }else{
    376465    // merge collapsed clusters
    377     while(merge(y,beta)){update(beta, tks,y);}
     466    while(merge(y,beta)){update1(beta, tks,y);}
    378467    if(fVerbose ){ cout << "dump after 1st merging " << endl;  dump(beta,y,tks);}
    379468  }
    … …  
    381470  // switch on outlier rejection
    382471  rho0=1./nt; for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){ k->pk =1.; }  // democratic
    383   niter=0; while((update(beta, tks,y,rho0) > 1.e-8)  && (niter++ < fMaxIterations)){  }
     472  niter=0; while((update2(beta, tks,y,rho0, fDzCutOff) > 1.e-8)  && (niter++ < fMaxIterations)){  }
    384473  if(fVerbose  ){ cout << "rho0=" << rho0 <<   " niter=" << niter <<  endl; dump(beta,y,tks);}
    385474
    … …  
    392481  // continue from freeze-out to Tstop (=1) without splitting, eliminate insignificant vertices
    393482  while(beta<=fBetaStop){
    394     while(purge(y,tks,rho0, beta)){
    395       niter=0; while((update(beta, tks, y, rho0) > 1.e-6)  && (niter++ < fMaxIterations)){  }
     483    while(purge(y,tks,rho0, beta, fDzCutOff)){
     484      niter=0; while((update2(beta, tks, y, rho0, fDzCutOff) > 1.e-6)  && (niter++ < fMaxIterations)){  }
    396485    }
    397486    beta/=fCoolingFactor;
    398     niter=0; while((update(beta, tks, y, rho0) > 1.e-6)  && (niter++ < fMaxIterations)){  }
     487    niter=0; while((update2(beta, tks, y, rho0, fDzCutOff) > 1.e-6)  && (niter++ < fMaxIterations)){  }
    399488  }
    400489
    … …  
    402491//   // new, one last round of cleaning at T=Tstop
    403492//   while(purge(y,tks,rho0, beta)){
    404 //     niter=0; while((update(beta, tks,y,rho0) > 1.e-6)  && (niter++ < fMaxIterations)){  }
     493//     niter=0; while((update2(beta, tks,y,rho0, fDzCutOff) > 1.e-6)  && (niter++ < fMaxIterations)){  }
    405494//   }
    406495
    … …  
    440529      const double invdt = 1.0/std::sqrt(tks[i].dt2);
    441530      if(tks[i].Z>0){
    442         double p = k->pk * exp(-beta*Eik(tks[i],*k)) / tks[i].Z;
    443         if( (tks[i].pi>0) && ( p > 0.5 ) ){
     531  double p = k->pk * exp(-beta*Eik(tks[i],*k)) / tks[i].Z;
     532  if( (tks[i].pi>0) && ( p > 0.5 ) ){
    444533          //std::cout << "pushing back " << i << ' ' << tks[i].tt << std::endl;
    445534          //vertexTracks.push_back(*(tks[i].tt)); tks[i].Z=0;
    … …  
    482571//------------------------------------------------------------------------------
    483572
    484 vector<VertexFinderDA4D::track_t>
    485 VertexFinderDA4D::fill()const{
    486 
    487   Candidate *candidate;
    488   track_t tr;
    489   Double_t z, dz, t, l, dt, d0, d0error;
    490 
    491   // prepare track data for clustering
    492   vector< track_t > tks;
    493 
    494   // loop over input tracks
    495   fItInputArray->Reset();
    496   while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    497   {
    498      //TBC everything in cm
    499      z = candidate->DZ/10;
    500      tr.z = z;
    501      dz = candidate->ErrorDZ/10;
    502      tr.dz2 = dz*dz          // track error
    503    // TBC: beamspot size induced error, take 0 for now.
    504    //   + (std::pow(beamspot.BeamWidthX()*cos(phi),2.)+std::pow(beamspot.BeamWidthY()*sin(phi),2.))/std::pow(tantheta,2.) // beam-width induced
    505       + fVertexSpaceSize*fVertexSpaceSize;                     // intrinsic vertex size, safer for outliers and short lived decays
    506 
    507      // TBC: the time is in ns for now TBC
    508      //t = candidate->Position.T()/c_light;
    509      t = candidate->InitialPosition.T()/c_light;
    510      l = candidate->L/c_light;
    511 
    512      double pt = candidate->Momentum.Pt();
    513      double eta = candidate->Momentum.Eta();
    514      double phi = candidate->Momentum.Phi();
    515 
    516      tr.pt = pt;
    517      tr.eta = eta;
    518      tr.phi = phi;
    519      tr.t = t; //
    520      tr.dtz = 0.;
    521      dt = candidate->ErrorT/c_light;
    522      tr.dt2 = dt*dt + fVertexTimeSize*fVertexTimeSize;   // the ~injected~ timing error plus a small minimum vertex size in time
    523      if (fD0CutOff>0){
    524 
    525      d0 = TMath::Abs(candidate->D0)/10.0;
    526      d0error = candidate->ErrorD0/10.0;
    527 
    528      tr.pi=1./(1.+exp((d0*d0)/(d0error*d0error) - fD0CutOff*fD0CutOff));  // reduce weight for high ip tracks
    529 
    530     }
    531     else{
    532       tr.pi=1.;
    533     }
    534     tr.tt=&(*candidate);
    535     tr.Z=1.;
    536 
    537     // TBC now putting track selection here (> fPTMin)
    538     if( tr.pi > 1e-3 && tr.pt > fMinPT) {
    539       tks.push_back(tr);
    540     }
    541 
    542   }
    543 
    544   return tks;
    545 }
    546 
    547 
    548 
    549 
    550 //------------------------------------------------------------------------------
    551 
    552 double VertexFinderDA4D::Eik(const track_t & t, const vertex_t &k ) const {
     573static double Eik(const track_t & t, const vertex_t &k)
     574{
    553575  return std::pow(t.z-k.z,2.)/t.dz2 + std::pow(t.t - k.t,2.)/t.dt2;
    554576}
    555577
    556 //--------------------------------------------------------------------------------
    557 
    558 void VertexFinderDA4D::dump(const double beta, const vector<vertex_t> & y, const vector<track_t> & tks0)const{
    559 
     578//------------------------------------------------------------------------------
     579
     580static void dump(const double beta, const vector<vertex_t> &y, const vector<track_t> &tks0)
     581{
    560582  // copy and sort for nicer printout
    561583  vector<track_t> tks;
    … …  
    564586
    565587  cout << "-----DAClusterizerInZT::dump ----" << endl;
    566   cout << " beta=" << beta << "   betamax= " << fBetaMax << endl;
     588  cout << " beta=" << beta << endl;
    567589  cout << "                                                               z= ";
    568590  cout.precision(4);
    … …  
    587609  cout  << endl;
    588610
    589   if(fVerbose){
    590     double E=0, F=0;
    591     cout << endl;
    592     cout << "----       z +/- dz        t +/- dt        ip +/-dip       pt    phi  eta    weights  ----" << endl;
    593     cout.precision(4);
    594     for(unsigned int i=0; i<tks.size(); i++){
    595       if (tks[i].Z>0){  F-=log(tks[i].Z)/beta;}
    596       double tz= tks[i].z;
    597       double tt= tks[i].t;
    598       //cout <<  setw (3)<< i << ")" <<  setw (8) << fixed << setprecision(4)<<  tz << " +/-" <<  setw (6)<< sqrt(tks[i].dz2)
    599       //     << setw(8) << fixed << setprecision(4) << tt << " +/-" << setw(6) << std::sqrt(tks[i].dt2)  ;
    600 
    601       double sump=0.;
    602       for(vector<vertex_t>::const_iterator k=y.begin(); k!=y.end(); k++){
    603         if((tks[i].pi>0)&&(tks[i].Z>0)){
    604           //double p=pik(beta,tks[i],*k);
    605           double p=k->pk * std::exp(-beta*Eik(tks[i],*k)) / tks[i].Z;
    606           if( p > 0.0001){
    607             //cout <<  setw (8) <<  setprecision(3) << p;
    608           }else{
    609             cout << "    .   ";
    610           }
    611           E+=p*Eik(tks[i],*k);
    612           sump+=p;
    613         }else{
    614             cout << "        ";
    615         }
     611  double E=0, F=0;
     612  cout << endl;
     613  cout << "----       z +/- dz        t +/- dt        ip +/-dip       pt    phi  eta    weights  ----" << endl;
     614  cout.precision(4);
     615  for(unsigned int i=0; i<tks.size(); i++){
     616    if (tks[i].Z>0){  F-=log(tks[i].Z)/beta;}
     617    double tz= tks[i].z;
     618    double tt= tks[i].t;
     619    //cout <<  setw (3)<< i << ")" <<  setw (8) << fixed << setprecision(4)<<  tz << " +/-" <<  setw (6)<< sqrt(tks[i].dz2)
     620    //     << setw(8) << fixed << setprecision(4) << tt << " +/-" << setw(6) << std::sqrt(tks[i].dt2)  ;
     621
     622    double sump=0.;
     623    for(vector<vertex_t>::const_iterator k=y.begin(); k!=y.end(); k++){
     624    if((tks[i].pi>0)&&(tks[i].Z>0)){
     625    //double p=pik(beta,tks[i],*k);
     626    double p=k->pk * std::exp(-beta*Eik(tks[i],*k)) / tks[i].Z;
     627    if( p > 0.0001){
     628      //cout <<  setw (8) <<  setprecision(3) << p;
     629    }else{
     630      cout << "    .   ";
     631    }
     632    E+=p*Eik(tks[i],*k);
     633    sump+=p;
     634  }else{
     635      cout << "        ";
     636  }
    616637      }
    617638      cout << endl;
    618639    }
    619640    cout << endl << "T=" << 1/beta  << " E=" << E << " n="<< y.size() << "  F= " << F <<  endl <<  "----------" << endl;
    620   }
    621 }
    622 
    623 
    624 
    625 
    626 
    627 //------------------------------------------------------------------------------
    628 
    629 
    630 double VertexFinderDA4D::update( double beta,
    631                                   vector<track_t> & tks,
    632                                   vector<vertex_t> & y ) const {
     641}
     642
     643//------------------------------------------------------------------------------
     644
     645static double update1(double beta, vector<track_t> &tks, vector<vertex_t> &y)
     646{
    633647  //update weights and vertex positions
    634648  // mass constrained annealing without noise
    … …  
    661675      // accumulate weighted z and weights for vertex update
    662676      for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); ++k){
    663         k->se  += tks[i].pi* k->ei / Zi;
    664         const double w = k->pk * tks[i].pi* k->ei / ( Zi * ( tks[i].dz2 * tks[i].dt2 ) );
    665         k->sw  += w;
    666         k->swz += w * tks[i].z;
     677  k->se  += tks[i].pi* k->ei / Zi;
     678  const double w = k->pk * tks[i].pi* k->ei / ( Zi * ( tks[i].dz2 * tks[i].dt2 ) );
     679  k->sw  += w;
     680  k->swz += w * tks[i].z;
    667681        k->swt += w * tks[i].t;
    668         k->swE += w * Eik(tks[i],*k);
     682  k->swE += w * Eik(tks[i],*k);
    669683      }
    670684    }else{
    … …  
    687701      k->Tc = 2.*k->swE/k->sw;
    688702    }else{
    689       if(fVerbose){cout << " a cluster melted away ?  pk=" << k->pk <<  " sumw=" << k->sw <<  endl;}
     703      // cout << " a cluster melted away ?  pk=" << k->pk <<  " sumw=" << k->sw <<  endl
    690704      k->Tc=-1;
    691705    }
    … …  
    700714//------------------------------------------------------------------------------
    701715
    702 double VertexFinderDA4D::update( double beta,
    703                                   vector<track_t> & tks,
    704                                   vector<vertex_t> & y,
    705                                   double & rho0 ) const {
     716static double update2(double beta, vector<track_t> &tks, vector<vertex_t> &y, double &rho0, double dzCutOff)
     717{
    706718  // MVF style, no more vertex weights, update tracks weights and vertex positions, with noise
    707719  // returns the squared sum of changes of vertex positions
    … …  
    720732
    721733    // update pik and Zi and Ti
    722     double Zi = rho0*std::exp(-beta*(fDzCutOff*fDzCutOff));// cut-off (eventually add finite size in time)
     734    double Zi = rho0*std::exp(-beta*(dzCutOff*dzCutOff));// cut-off (eventually add finite size in time)
    723735    //double Ti = 0.; // dt0*std::exp(-beta*fDtCutOff);
    724736    for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
    … …  
    732744       // accumulate weighted z and weights for vertex update
    733745      for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
    734         k->se += tks[i].pi* k->ei / Zi;
    735         double w = k->pk * tks[i].pi * k->ei /( Zi * ( tks[i].dz2 * tks[i].dt2 ) );
    736         k->sw  += w;
    737         k->swz += w * tks[i].z;
     746  k->se += tks[i].pi* k->ei / Zi;
     747  double w = k->pk * tks[i].pi * k->ei /( Zi * ( tks[i].dz2 * tks[i].dt2 ) );
     748  k->sw  += w;
     749  k->swz += w * tks[i].z;
    738750        k->swt += w * tks[i].t;
    739         k->swE += w * Eik(tks[i],*k);
     751  k->swE += w * Eik(tks[i],*k);
    740752      }
    741753    }
    … …  
    754766      k->Tc  = 2*k->swE/k->sw;
    755767    }else{
    756       if(fVerbose){cout << " a cluster melted away ?  pk=" << k->pk <<  " sumw=" << k->sw <<  endl;}
     768      // cout << " a cluster melted away ?  pk=" << k->pk <<  " sumw=" << k->sw <<  endl;
    757769      k->Tc = 0;
    758770    }
    … …  
    766778//------------------------------------------------------------------------------
    767779
    768 bool VertexFinderDA4D::merge(vector<vertex_t> & y)const{
     780static bool merge(vector<vertex_t> &y)
     781{
    769782  // merge clusters that collapsed or never separated, return true if vertices were merged, false otherwise
    770783
    … …  
    794807//------------------------------------------------------------------------------
    795808
    796 bool VertexFinderDA4D::merge(vector<vertex_t> & y, double & beta)const{
     809static bool merge(vector<vertex_t> &y, double &beta)
     810{
    797811  // merge clusters that collapsed or never separated,
    798812  // only merge if the estimated critical temperature of the merged vertex is below the current temperature
    … …  
    809823
    810824      if(Tc*beta<1){
    811         if(rho>0){
    812           k->z = ( k->pk * k->z + (k+1)->z * (k+1)->pk)/rho;
     825  if(rho>0){
     826    k->z = ( k->pk * k->z + (k+1)->z * (k+1)->pk)/rho;
    813827          k->t = ( k->pk * k->t + (k+1)->t * (k+1)->pk)/rho;
    814         }else{
    815           k->z = 0.5*(k->z + (k+1)->z);
     828  }else{
     829    k->z = 0.5*(k->z + (k+1)->z);
    816830          k->t = 0.5*(k->t + (k+1)->t);
    817         }
    818         k->pk  = rho;
    819         k->sw += (k+1)->sw;
    820         k->swE = swE;
    821         k->Tc  = Tc;
    822         y.erase(k+1);
    823         return true;
     831  }
     832  k->pk  = rho;
     833  k->sw += (k+1)->sw;
     834  k->swE = swE;
     835  k->Tc  = Tc;
     836  y.erase(k+1);
     837  return true;
    824838      }
    825839    }
    … …  
    831845//------------------------------------------------------------------------------
    832846
    833 bool VertexFinderDA4D::purge(vector<vertex_t> & y, vector<track_t> & tks, double & rho0, const double beta)const{
     847static bool purge(vector<vertex_t> &y, vector<track_t> &tks, double & rho0, const double beta, const double dzCutOff)
     848{
    834849  // eliminate clusters with only one significant/unique track
    835850  if(y.size()<2)  return false;
    … …  
    841856    int nUnique=0;
    842857    double sump=0;
    843     double pmax=k->pk/(k->pk+rho0*exp(-beta*fDzCutOff*fDzCutOff));
     858    double pmax=k->pk/(k->pk+rho0*exp(-beta*dzCutOff*dzCutOff));
    844859    for(unsigned int i=0; i<nt; i++){
    845860      if(tks[i].Z > 0){
    846         double p = k->pk * std::exp(-beta*Eik(tks[i],*k)) / tks[i].Z ;
    847         sump+=p;
    848         if( (p > 0.9*pmax) && (tks[i].pi>0) ){ nUnique++; }
     861  double p = k->pk * std::exp(-beta*Eik(tks[i],*k)) / tks[i].Z ;
     862  sump+=p;
     863  if( (p > 0.9*pmax) && (tks[i].pi>0) ){ nUnique++; }
    849864      }
    850865    }
    … …  
    857872
    858873  if(k0!=y.end()){
    859     if(fVerbose){cout << "eliminating prototype at " << k0->z << "," << k0->t
    860                       << " with sump=" << sumpmin << endl;}
     874    //cout << "eliminating prototype at " << k0->z << "," << k0->t << " with sump=" << sumpmin << endl;
    861875    //rho0+=k0->pk;
    862876    y.erase(k0);
    … …  
    869883//------------------------------------------------------------------------------
    870884
    871 double VertexFinderDA4D::beta0( double betamax,
    872                                  vector<track_t> & tks,
    873                                  vector<vertex_t> & y ) const {
     885static double beta0(double betamax, vector<track_t> &tks, vector<vertex_t> &y, const double coolingFactor)
     886{
    874887
    875888  double T0=0;  // max Tc for beta=0
    … …  
    906919
    907920  if (T0>1./betamax){
    908     return betamax/pow(fCoolingFactor, int(std::log(T0*betamax)/std::log(fCoolingFactor))-1 );
     921    return betamax/pow(coolingFactor, int(std::log(T0*betamax)/std::log(coolingFactor))-1 );
    909922  }else{
    910923    // ensure at least one annealing step
    911     return betamax/fCoolingFactor;
    912   }
    913 }
    914 
    915 //------------------------------------------------------------------------------
    916 
    917 bool VertexFinderDA4D::split( double beta,
    918                                vector<track_t> & tks,
    919                                vector<vertex_t> & y) const {
     924    return betamax/coolingFactor;
     925  }
     926}
     927
     928//------------------------------------------------------------------------------
     929
     930static bool split(double beta, vector<track_t> &tks, vector<vertex_t> &y)
     931{
    920932  // split only critical vertices (Tc >~ T=1/beta   <==>   beta*Tc>~1)
    921933  // an update must have been made just before doing this (same beta, no merging)
    922934  // returns true if at least one cluster was split
    923935
    924   constexpr double epsilon=1e-3;      // split all single vertices by 10 um
    925   bool split=false;
     936  const double epsilon = 1e-3; // split all single vertices by 10 um
     937  bool split = false;
    926938
    927939  // avoid left-right biases by splitting highest Tc first
    … …  
    943955    for(unsigned int i=0; i<tks.size(); i++){
    944956      if(tks[i].Z>0){
    945         //sumpi+=tks[i].pi;
    946         double p=y[ik].pk * exp(-beta*Eik(tks[i],y[ik])) / tks[i].Z*tks[i].pi;
    947         double w=p/(tks[i].dz2 * tks[i].dt2);
    948         if(tks[i].z < y[ik].z){
    949           p1+=p; z1+=w*tks[i].z; t1+=w*tks[i].t; w1+=w;
    950         }else{
    951           p2+=p; z2+=w*tks[i].z; t2+=w*tks[i].t; w2+=w;
    952         }
     957  //sumpi+=tks[i].pi;
     958  double p=y[ik].pk * exp(-beta*Eik(tks[i],y[ik])) / tks[i].Z*tks[i].pi;
     959  double w=p/(tks[i].dz2 * tks[i].dt2);
     960  if(tks[i].z < y[ik].z){
     961    p1+=p; z1+=w*tks[i].z; t1+=w*tks[i].t; w1+=w;
     962  }else{
     963    p2+=p; z2+=w*tks[i].z; t2+=w*tks[i].t; w2+=w;
     964  }
    953965      }
    954966    }
    … …  
    974986     // adjust remaining pointers
    975987      for(unsigned int jc=ic; jc<critical.size(); jc++){
    976         if (critical[jc].second>ik) {critical[jc].second++;}
     988        if (critical[jc].second>ik) {critical[jc].second++;}
    977989      }
    978990    }
    … …  
    985997//------------------------------------------------------------------------------
    986998
    987 void VertexFinderDA4D::splitAll( vector<vertex_t> & y ) const {
    988 
    989 
    990   constexpr double epsilon=1e-3;      // split all single vertices by 10 um
    991   constexpr double zsep=2*epsilon;    // split vertices that are isolated by at least zsep (vertices that haven't collapsed)
    992   constexpr double tsep=2*epsilon;    // check t as well
     999void splitAll(vector<vertex_t> &y)
     1000{
     1001
     1002
     1003  const double epsilon=1e-3;      // split all single vertices by 10 um
     1004  const double zsep=2*epsilon;    // split vertices that are isolated by at least zsep (vertices that haven't collapsed)
     1005  const double tsep=2*epsilon;    // check t as well
    9931006
    9941007  vector<vertex_t> y1;

  • modules/VertexFinderDA4D.h

    r95b4e9f r4154bbd  
    3030  void Finish();
    3131
    32   struct track_t
    33   {
    34     double z;      // z-coordinate at point of closest approach to the beamline
    35     double t;      // t-coordinate at point of closest approach to the beamline
    36     double dz2;    // square of the error of z(pca)
    37     double dtz;    // covariance of z-t
    38     double dt2;    // square of the error of t(pca)
    39     Candidate* tt; // a pointer to the Candidate Track
    40     double Z;      // Z[i]   for DA clustering
    41     double pi;     // track weight
    42     double pt;
    43     double eta;
    44     double phi;
    45   };
    46 
    47 
    48   struct vertex_t
    49   {
    50     double z;
    51     double t;
    52     double pk;     // vertex weight for "constrained" clustering
    53     // --- temporary numbers, used during update
    54     double ei;
    55     double sw;
    56     double swz;
    57     double swt;
    58     double se;
    59     // ---for Tc
    60     double swE;
    61     double Tc;
    62   };
    63 
    64   void clusterize(const TObjArray & tracks, TObjArray & clusters);
    65 
     32  void clusterize(const TObjArray &tracks, TObjArray &clusters);
    6633  std::vector< Candidate* > vertices();
    67 
    68   std::vector<track_t> fill() const;
    69 
    70   bool split(double beta,
    71              std::vector<track_t> & tks,
    72              std::vector<vertex_t> & y) const;
    73 
    74   double update(double beta,
    75                 std::vector<track_t> & tks,
    76                 std::vector<vertex_t> & y) const;
    77 
    78   double update(double beta,
    79                 std::vector<track_t> & tks,
    80                 std::vector<vertex_t> & y,
    81                 double &)const;
    82 
    83   void dump(const double beta, const std::vector<vertex_t> & y, const std::vector<track_t> & tks) const;
    84   bool merge(std::vector<vertex_t> &) const;
    85   bool merge(std::vector<vertex_t> &, double &) const;
    86   bool purge(std::vector<vertex_t> &, std::vector<track_t> & , double &, const double) const;
    87 
    88   void splitAll(std::vector<vertex_t> &y) const;
    89 
    90   double beta0(const double betamax,
    91                std::vector<track_t> &tks,
    92                std::vector<vertex_t> &y)const;
    93 
    94   double Eik(const track_t &t, const vertex_t &k)const;
    9534
    9635private:

  • modules/VertexSorter.cc

    r95b4e9f r4154bbd  
    115115  Candidate *candidate, *jetCandidate, *beamSpotCandidate;
    116116  map<Int_t, UInt_t> clusterIDToIndex;
     117  map<Int_t, UInt_t>::const_iterator itClusterIDToIndex;
    117118  map<Int_t, Double_t> clusterIDToSumPT2;
     119  map<Int_t, Double_t>::const_iterator itClusterIDToSumPT2;
    118120  vector<pair<Int_t, Double_t> > sortedClusterIDs;
     121  vector<pair<Int_t, Double_t> >::const_iterator itSortedClusterIDs;
    119122
    120123  for (Int_t iCluster = 0; iCluster < fInputArray->GetEntries (); iCluster++)
    … …  
    161164        }
    162165
    163       for (const auto &clusterID : clusterIDToSumPT2)
    164         sortedClusterIDs.push_back (make_pair (clusterID.first, clusterID.second));
     166      for (itClusterIDToSumPT2 = clusterIDToSumPT2.begin(); itClusterIDToSumPT2 != clusterIDToSumPT2.end(); ++itClusterIDToSumPT2)
     167        sortedClusterIDs.push_back (make_pair (itClusterIDToSumPT2->first, itClusterIDToSumPT2->second));
    165168      sort (sortedClusterIDs.begin (), sortedClusterIDs.end (), secondDescending);
    166169    }
    … …  
    193196          if (candidate->IsPU)
    194197            continue;
    195           for (const auto &clusterID : clusterIDToIndex)
     198          for (itClusterIDToIndex = clusterIDToIndex.begin(); itClusterIDToIndex != clusterIDToIndex.end(); ++itClusterIDToIndex)
    196199            {
    197               if (candidate->ClusterIndex != clusterID.first)
     200              if (candidate->ClusterIndex != itClusterIDToIndex->first)
    198201                continue;
    199               clusterIDToSumPT2.at (clusterID.first) += candidate->Momentum.Pt () * candidate->Momentum.Pt ();
     202              clusterIDToSumPT2.at (itClusterIDToIndex->first) += candidate->Momentum.Pt () * candidate->Momentum.Pt ();
    200203            }
    201204        }
    202205
    203       for (const auto &clusterID : clusterIDToSumPT2)
    204         sortedClusterIDs.push_back (make_pair (clusterID.first, clusterID.second));
     206      for (itClusterIDToSumPT2 = clusterIDToSumPT2.begin(); itClusterIDToSumPT2 != clusterIDToSumPT2.end(); ++itClusterIDToSumPT2)
     207        sortedClusterIDs.push_back (make_pair (itClusterIDToSumPT2->first, itClusterIDToSumPT2->second));
    205208      sort (sortedClusterIDs.begin (), sortedClusterIDs.end (), secondDescending);
    206209    }
    … …  
    214217      throw 0;
    215218    }
    216 
    217   for (const auto &clusterID : sortedClusterIDs)
    218     {
    219       Candidate *cluster = (Candidate *) fInputArray->At (clusterIDToIndex.at (clusterID.first));
     219  for (itSortedClusterIDs = sortedClusterIDs.begin(); itSortedClusterIDs != sortedClusterIDs.end(); ++itSortedClusterIDs)
     220    {
     221      Candidate *cluster = (Candidate *) fInputArray->At (clusterIDToIndex.at (itSortedClusterIDs->first));
    220222      if (fMethod == "BTV")
    221         cluster->BTVSumPT2 = clusterID.second;
     223        cluster->BTVSumPT2 = itSortedClusterIDs->second;
    222224      else if (fMethod == "GenClosest")
    223         cluster->GenDeltaZ = clusterID.second;
     225        cluster->GenDeltaZ = itSortedClusterIDs->second;
    224226      else if (fMethod == "GenBest")
    225         cluster->GenSumPT2 = clusterID.second;
     227        cluster->GenSumPT2 = itSortedClusterIDs->second;
    226228      fOutputArray->Add (cluster);
    227229    }
Note: See TracChangeset for help on using the changeset viewer.