Fork me on GitHub

source: git/modules/VertexFinderDA4D.cc @ 77e9ae1

ImprovedOutputFileTimingllp
Last change on this file since 77e9ae1 was 77e9ae1, checked in by Pavel Demin <pavel-demin@…>, 18 months ago

set Standard to Cpp03 in .clang-format

  • Property mode set to 100644
File size: 33.7 KB
Line 
1/** \class VertexFinderDA4D
2 *
3 *  Cluster vertices from tracks using deterministic annealing and timing information
4 *
5 *  \authors M. Selvaggi, L. Gray
6 *
7 */
8
9#include "modules/VertexFinderDA4D.h"
10#include "classes/DelphesClasses.h"
11#include "classes/DelphesFactory.h"
12#include "classes/DelphesFormula.h"
13#include "classes/DelphesPileUpReader.h"
14
15#include "ExRootAnalysis/ExRootClassifier.h"
16#include "ExRootAnalysis/ExRootFilter.h"
17#include "ExRootAnalysis/ExRootResult.h"
18
19#include "TDatabasePDG.h"
20#include "TFormula.h"
21#include "TLorentzVector.h"
22#include "TMath.h"
23#include "TMatrixT.h"
24#include "TObjArray.h"
25#include "TRandom3.h"
26#include "TString.h"
27#include "TVector3.h"
28
29#include <algorithm>
30#include <iostream>
31#include <stdexcept>
32#include <utility>
33#include <vector>
34
35using namespace std;
36
37static const Double_t mm = 1.;
38static const Double_t m = 1000. * mm;
39static const Double_t ns = 1.;
40static const Double_t s = 1.e+9 * ns;
41static const Double_t c_light = 2.99792458e+8 * m / s;
42
43struct track_t
44{
45  double z; // z-coordinate at point of closest approach to the beamline
46  double t; // t-coordinate at point of closest approach to the beamline
47  double dz2; // square of the error of z(pca)
48  double dtz; // covariance of z-t
49  double dt2; // square of the error of t(pca)
50  Candidate *tt; // a pointer to the Candidate Track
51  double Z; // Z[i]   for DA clustering
52  double pi; // track weight
53  double pt;
54  double eta;
55  double phi;
56};
57
58struct vertex_t
59{
60  double z;
61  double t;
62  double pk; // vertex weight for "constrained" clustering
63  // --- temporary numbers, used during update
64  double ei;
65  double sw;
66  double swz;
67  double swt;
68  double se;
69  // ---for Tc
70  double swE;
71  double Tc;
72};
73
74static bool split(double beta, std::vector<track_t> &tks, std::vector<vertex_t> &y);
75static double update1(double beta, std::vector<track_t> &tks, std::vector<vertex_t> &y);
76static double update2(double beta, std::vector<track_t> &tks, std::vector<vertex_t> &y, double &rho0, const double dzCutOff);
77static void dump(const double beta, const std::vector<vertex_t> &y, const std::vector<track_t> &tks);
78static bool merge(std::vector<vertex_t> &);
79static bool merge(std::vector<vertex_t> &, double &);
80static bool purge(std::vector<vertex_t> &, std::vector<track_t> &, double &, const double, const double);
81static void splitAll(std::vector<vertex_t> &y);
82static double beta0(const double betamax, std::vector<track_t> &tks, std::vector<vertex_t> &y, const double coolingFactor);
83static double Eik(const track_t &t, const vertex_t &k);
84
85static bool recTrackLessZ1(const track_t &tk1, const track_t &tk2)
86{
87  return tk1.z < tk2.z;
88}
89
90using namespace std;
91
92//------------------------------------------------------------------------------
93
94VertexFinderDA4D::VertexFinderDA4D() :
95  fVerbose(0), fMinPT(0), fVertexSpaceSize(0), fVertexTimeSize(0),
96  fUseTc(0), fBetaMax(0), fBetaStop(0), fCoolingFactor(0),
97  fMaxIterations(0), fDzCutOff(0), fD0CutOff(0), fDtCutOff(0)
98{
99}
100
101//------------------------------------------------------------------------------
102
103VertexFinderDA4D::~VertexFinderDA4D()
104{
105}
106
107//------------------------------------------------------------------------------
108
109void VertexFinderDA4D::Init()
110{
111
112  fVerbose = GetBool("Verbose", 1);
113  fMinPT = GetDouble("MinPT", 0.1);
114  fVertexSpaceSize = GetDouble("VertexSpaceSize", 0.5); //in mm
115  fVertexTimeSize = GetDouble("VertexTimeSize", 10E-12); //in s
116  fUseTc = GetBool("UseTc", 1);
117  fBetaMax = GetDouble("BetaMax ", 0.1);
118  fBetaStop = GetDouble("BetaStop", 1.0);
119  fCoolingFactor = GetDouble("CoolingFactor", 0.8);
120  fMaxIterations = GetInt("MaxIterations", 100);
121  fDzCutOff = GetDouble("DzCutOff", 40); // Adaptive Fitter uses 30 mm but that appears to be a bit tight here sometimes
122  fD0CutOff = GetDouble("D0CutOff", 30);
123  fDtCutOff = GetDouble("DtCutOff", 100E-12); // dummy
124
125  // convert stuff in cm, ns
126  fVertexSpaceSize /= 10.0;
127  fVertexTimeSize *= 1E9;
128  fDzCutOff /= 10.0; // Adaptive Fitter uses 3.0 but that appears to be a bit tight here sometimes
129  fD0CutOff /= 10.0;
130
131  fInputArray = ImportArray(GetString("InputArray", "TrackSmearing/tracks"));
132  fItInputArray = fInputArray->MakeIterator();
133
134  fOutputArray = ExportArray(GetString("OutputArray", "tracks"));
135  fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
136}
137
138//------------------------------------------------------------------------------
139
140void VertexFinderDA4D::Finish()
141{
142  if(fItInputArray) delete fItInputArray;
143}
144
145//------------------------------------------------------------------------------
146
147void VertexFinderDA4D::Process()
148{
149  Candidate *candidate, *track;
150  TObjArray *ClusterArray;
151  ClusterArray = new TObjArray;
152  TIterator *ItClusterArray;
153  Int_t ivtx = 0;
154
155  fInputArray->Sort();
156
157  TLorentzVector pos, mom;
158  if(fVerbose)
159  {
160    cout << " start processing vertices ..." << endl;
161    cout << " Found " << fInputArray->GetEntriesFast() << " input tracks" << endl;
162    //loop over input tracks
163    fItInputArray->Reset();
164    while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
165    {
166      pos = candidate->InitialPosition;
167      mom = candidate->Momentum;
168
169      cout << "pt: " << mom.Pt() << ", eta: " << mom.Eta() << ", phi: " << mom.Phi() << ", z: " << candidate->DZ / 10 << endl;
170    }
171  }
172
173  // clusterize tracks in Z
174  clusterize(*fInputArray, *ClusterArray);
175
176  if(fVerbose)
177  {
178    std::cout << " clustering returned  " << ClusterArray->GetEntriesFast() << " clusters  from " << fInputArray->GetEntriesFast() << " selected tracks" << std::endl;
179  }
180
181  //loop over vertex candidates
182  ItClusterArray = ClusterArray->MakeIterator();
183  ItClusterArray->Reset();
184  while((candidate = static_cast<Candidate *>(ItClusterArray->Next())))
185  {
186
187    double meantime = 0.;
188    double expv_x2 = 0.;
189    double normw = 0.;
190    double errtime = 0;
191
192    double meanpos = 0.;
193    double meanerr2 = 0.;
194    double normpos = 0.;
195    double errpos = 0.;
196
197    double sumpt2 = 0.;
198
199    int itr = 0;
200
201    if(fVerbose) cout << "this vertex has: " << candidate->GetCandidates()->GetEntriesFast() << " tracks" << endl;
202
203    // loop over tracks belonging to this vertex
204    TIter it1(candidate->GetCandidates());
205    it1.Reset();
206
207    while((track = static_cast<Candidate *>(it1.Next())))
208    {
209
210      itr++;
211      // TBC: the time is in ns for now TBC
212      double t = track->InitialPosition.T() / c_light;
213      double dt = track->ErrorT / c_light;
214      const double time = t;
215      const double inverr = 1.0 / dt;
216      meantime += time * inverr;
217      expv_x2 += time * time * inverr;
218      normw += inverr;
219
220      // compute error position TBC
221      const double pt = track->Momentum.Pt();
222      const double z = track->DZ / 10.0;
223      const double err_pt = track->ErrorPT;
224      const double err_z = track->ErrorDZ;
225
226      const double wi = (pt / (err_pt * err_z)) * (pt / (err_pt * err_z));
227      meanpos += z * wi;
228
229      meanerr2 += err_z * err_z * wi;
230      normpos += wi;
231      sumpt2 += pt * pt;
232
233      // while we are here store cluster index in tracks
234      track->ClusterIndex = ivtx;
235    }
236
237    meantime = meantime / normw;
238    expv_x2 = expv_x2 / normw;
239    errtime = TMath::Sqrt((expv_x2 - meantime * meantime) / itr);
240    meanpos = meanpos / normpos;
241    meanerr2 = meanerr2 / normpos;
242    errpos = TMath::Sqrt(meanerr2 / itr);
243
244    candidate->Position.SetXYZT(0.0, 0.0, meanpos * 10.0, meantime * c_light);
245    candidate->PositionError.SetXYZT(0.0, 0.0, errpos * 10.0, errtime * c_light);
246    candidate->SumPT2 = sumpt2;
247    candidate->ClusterNDF = itr;
248    candidate->ClusterIndex = ivtx;
249
250    fVertexOutputArray->Add(candidate);
251
252    ivtx++;
253
254    if(fVerbose)
255    {
256      std::cout << "x,y,z";
257      std::cout << ",t";
258      std::cout << "=" << candidate->Position.X() / 10.0 << " " << candidate->Position.Y() / 10.0 << " " << candidate->Position.Z() / 10.0;
259      std::cout << " " << candidate->Position.T() / c_light;
260
261      std::cout << std::endl;
262      std::cout << "sumpt2 " << candidate->SumPT2 << endl;
263
264      std::cout << "ex,ey,ez";
265      std::cout << ",et";
266      std::cout << "=" << candidate->PositionError.X() / 10.0 << " " << candidate->PositionError.Y() / 10.0 << " " << candidate->PositionError.Z() / 10.0;
267      std::cout << " " << candidate->PositionError.T() / c_light;
268      std::cout << std::endl;
269    }
270  } // end of cluster loop
271
272  if(fVerbose)
273  {
274    std::cout << "PrimaryVertexProducerAlgorithm::vertices candidates =" << ClusterArray->GetEntriesFast() << std::endl;
275  }
276
277  //TBC maybe this can be done later
278  // sort vertices by pt**2  vertex (aka signal vertex tagging)
279  /*if(pvs.size()>1){
280      sort(pvs.begin(), pvs.end(), VertexHigherPtSquared());
281    }
282     */
283
284  delete ClusterArray;
285}
286
287//------------------------------------------------------------------------------
288
289void VertexFinderDA4D::clusterize(const TObjArray &tracks, TObjArray &clusters)
290{
291  if(fVerbose)
292  {
293    cout << "###################################################" << endl;
294    cout << "# VertexFinderDA4D::clusterize   nt=" << tracks.GetEntriesFast() << endl;
295    cout << "###################################################" << endl;
296  }
297
298  vector<Candidate *> pv = vertices();
299
300  if(fVerbose)
301  {
302    cout << "# VertexFinderDA4D::clusterize   pv.size=" << pv.size() << endl;
303  }
304  if(pv.size() == 0)
305  {
306    return;
307  }
308
309  // convert into vector of candidates
310  //TObjArray *ClusterArray = pv.begin()->GetCandidates();
311  //Candidate *aCluster = static_cast<Candidate*>(&(pv.at(0)));
312  Candidate *aCluster = pv.at(0);
313
314  // fill into clusters and merge
315
316  if(fVerbose)
317  {
318    std::cout << '\t' << 0;
319    std::cout << ' ' << (*pv.begin())->Position.Z() / 10.0 << ' ' << (*pv.begin())->Position.T() / c_light << std::endl;
320  }
321
322  for(vector<Candidate *>::iterator k = pv.begin() + 1; k != pv.end(); k++)
323  {
324    if(fVerbose)
325    {
326      std::cout << '\t' << std::distance(pv.begin(), k);
327      std::cout << ' ' << (*k)->Position.Z() << ' ' << (*k)->Position.T() << std::endl;
328    }
329
330    // TBC - check units here
331    if(std::abs((*k)->Position.Z() - (*(k - 1))->Position.Z()) / 10.0 > (2 * fVertexSpaceSize) || std::abs((*k)->Position.T() - (*(k - 1))->Position.Z()) / c_light > 2 * 0.010)
332    {
333      // close a cluster
334      clusters.Add(aCluster);
335      //aCluster.clear();
336    }
337    //for(unsigned int i=0; i<k->GetCandidates().GetEntriesFast(); i++){
338    aCluster = *k;
339    //}
340  }
341  clusters.Add(aCluster);
342
343  if(fVerbose)
344  {
345    std::cout << "# VertexFinderDA4D::clusterize clusters.size=" << clusters.GetEntriesFast() << std::endl;
346  }
347}
348
349//------------------------------------------------------------------------------
350
351vector<Candidate *> VertexFinderDA4D::vertices()
352{
353  Candidate *candidate;
354  UInt_t clusterIndex = 0;
355  vector<Candidate *> clusters;
356
357  vector<track_t> tks;
358  track_t tr;
359  Double_t z, dz, t, l, dt, d0, d0error;
360
361  // loop over input tracks
362  fItInputArray->Reset();
363  while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
364  {
365    //TBC everything in cm
366    z = candidate->DZ / 10;
367    tr.z = z;
368    dz = candidate->ErrorDZ / 10;
369    tr.dz2 = dz * dz // track error
370      //TBC: beamspot size induced error, take 0 for now.
371      // + (std::pow(beamspot.BeamWidthX()*cos(phi),2.)+std::pow(beamspot.BeamWidthY()*sin(phi),2.))/std::pow(tantheta,2.) // beam-width induced
372      + fVertexSpaceSize * fVertexSpaceSize; // intrinsic vertex size, safer for outliers and short lived decays
373
374    // TBC: the time is in ns for now TBC
375    //t = candidate->Position.T()/c_light;
376    t = candidate->InitialPosition.T() / c_light;
377    l = candidate->L / c_light;
378    double pt = candidate->Momentum.Pt();
379    double eta = candidate->Momentum.Eta();
380    double phi = candidate->Momentum.Phi();
381
382    tr.pt = pt;
383    tr.eta = eta;
384    tr.phi = phi;
385    tr.t = t; //
386    tr.dtz = 0.;
387    dt = candidate->ErrorT / c_light;
388    tr.dt2 = dt * dt + fVertexTimeSize * fVertexTimeSize; // the ~injected~ timing error plus a small minimum vertex size in time
389    if(fD0CutOff > 0)
390    {
391
392      d0 = TMath::Abs(candidate->D0) / 10.0;
393      d0error = candidate->ErrorD0 / 10.0;
394
395      tr.pi = 1. / (1. + exp((d0 * d0) / (d0error * d0error) - fD0CutOff * fD0CutOff)); // reduce weight for high ip tracks
396    }
397    else
398    {
399      tr.pi = 1.;
400    }
401    tr.tt = &(*candidate);
402    tr.Z = 1.;
403
404    // TBC now putting track selection here (> fPTMin)
405    if(tr.pi > 1e-3 && tr.pt > fMinPT)
406    {
407      tks.push_back(tr);
408    }
409  }
410
411  //print out input tracks
412
413  if(fVerbose)
414  {
415    std::cout << " start processing vertices ..." << std::endl;
416    std::cout << " Found " << tks.size() << " input tracks" << std::endl;
417    //loop over input tracks
418
419    for(std::vector<track_t>::const_iterator it = tks.begin(); it != tks.end(); it++)
420    {
421      double z = it->z;
422      double pt = it->pt;
423      double eta = it->eta;
424      double phi = it->phi;
425      double t = it->t;
426
427      std::cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", z: " << z << ", t: " << t << std::endl;
428    }
429  }
430
431  unsigned int nt = tks.size();
432  double rho0 = 0.0; // start with no outlier rejection
433
434  if(tks.empty()) return clusters;
435
436  vector<vertex_t> y; // the vertex prototypes
437
438  // initialize:single vertex at infinite temperature
439  vertex_t vstart;
440  vstart.z = 0.;
441  vstart.t = 0.;
442  vstart.pk = 1.;
443  y.push_back(vstart);
444  int niter = 0; // number of iterations
445
446  // estimate first critical temperature
447  double beta = beta0(fBetaMax, tks, y, fCoolingFactor);
448  niter = 0;
449  while((update1(beta, tks, y) > 1.e-6) && (niter++ < fMaxIterations))
450  {
451  }
452
453  // annealing loop, stop when T<Tmin  (i.e. beta>1/Tmin)
454  while(beta < fBetaMax)
455  {
456
457    if(fUseTc)
458    {
459      update1(beta, tks, y);
460      while(merge(y, beta))
461      {
462        update1(beta, tks, y);
463      }
464      split(beta, tks, y);
465      beta = beta / fCoolingFactor;
466    }
467    else
468    {
469      beta = beta / fCoolingFactor;
470      splitAll(y);
471    }
472
473    // make sure we are not too far from equilibrium before cooling further
474    niter = 0;
475    while((update1(beta, tks, y) > 1.e-6) && (niter++ < fMaxIterations))
476    {
477    }
478  }
479
480  if(fUseTc)
481  {
482    // last round of splitting, make sure no critical clusters are left
483    update1(beta, tks, y);
484    while(merge(y, beta))
485    {
486      update1(beta, tks, y);
487    }
488    unsigned int ntry = 0;
489    while(split(beta, tks, y) && (ntry++ < 10))
490    {
491      niter = 0;
492      while((update1(beta, tks, y) > 1.e-6) && (niter++ < fMaxIterations))
493      {
494      }
495      merge(y, beta);
496      update1(beta, tks, y);
497    }
498  }
499  else
500  {
501    // merge collapsed clusters
502    while(merge(y, beta))
503    {
504      update1(beta, tks, y);
505    }
506    if(fVerbose)
507    {
508      cout << "dump after 1st merging " << endl;
509      dump(beta, y, tks);
510    }
511  }
512
513  // switch on outlier rejection
514  rho0 = 1. / nt;
515  for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
516  {
517    k->pk = 1.;
518  } // democratic
519  niter = 0;
520  while((update2(beta, tks, y, rho0, fDzCutOff) > 1.e-8) && (niter++ < fMaxIterations))
521  {
522  }
523  if(fVerbose)
524  {
525    cout << "rho0=" << rho0 << " niter=" << niter << endl;
526    dump(beta, y, tks);
527  }
528
529  // merge again  (some cluster split by outliers collapse here)
530  while(merge(y))
531  {
532  }
533  if(fVerbose)
534  {
535    cout << "dump after 2nd merging " << endl;
536    dump(beta, y, tks);
537  }
538
539  // continue from freeze-out to Tstop (=1) without splitting, eliminate insignificant vertices
540  while(beta <= fBetaStop)
541  {
542    while(purge(y, tks, rho0, beta, fDzCutOff))
543    {
544      niter = 0;
545      while((update2(beta, tks, y, rho0, fDzCutOff) > 1.e-6) && (niter++ < fMaxIterations))
546      {
547      }
548    }
549    beta /= fCoolingFactor;
550    niter = 0;
551    while((update2(beta, tks, y, rho0, fDzCutOff) > 1.e-6) && (niter++ < fMaxIterations))
552    {
553    }
554  }
555
556  //   // new, one last round of cleaning at T=Tstop
557  //   while(purge(y,tks,rho0, beta)){
558  //     niter=0; while((update2(beta, tks,y,rho0, fDzCutOff) > 1.e-6)  && (niter++ < fMaxIterations)){  }
559  //   }
560
561  if(fVerbose)
562  {
563    cout << "Final result, rho0=" << rho0 << endl;
564    dump(beta, y, tks);
565  }
566
567  // select significant tracks and use a TransientVertex as a container
568  //GlobalError dummyError;
569
570  // ensure correct normalization of probabilities, should make double assginment reasonably impossible
571  for(unsigned int i = 0; i < nt; i++)
572  {
573    tks[i].Z = rho0 * exp(-beta * (fDzCutOff * fDzCutOff));
574    for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
575    {
576      tks[i].Z += k->pk * exp(-beta * Eik(tks[i], *k));
577    }
578  }
579
580  for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
581  {
582
583    DelphesFactory *factory = GetFactory();
584    candidate = factory->NewCandidate();
585
586    //cout<<"new vertex"<<endl;
587    //GlobalPoint pos(0, 0, k->z);
588    double time = k->t;
589    double z = k->z;
590    //vector< reco::TransientTrack > vertexTracks;
591    //double max_track_time_err2 = 0;
592    double mean = 0.;
593    double expv_x2 = 0.;
594    double normw = 0.;
595    for(unsigned int i = 0; i < nt; i++)
596    {
597      const double invdt = 1.0 / std::sqrt(tks[i].dt2);
598      if(tks[i].Z > 0)
599      {
600        double p = k->pk * exp(-beta * Eik(tks[i], *k)) / tks[i].Z;
601        if((tks[i].pi > 0) && (p > 0.5))
602        {
603          //std::cout << "pushing back " << i << ' ' << tks[i].tt << std::endl;
604          //vertexTracks.push_back(*(tks[i].tt)); tks[i].Z=0;
605
606          candidate->AddCandidate(tks[i].tt);
607          tks[i].Z = 0;
608
609          mean += tks[i].t * invdt * p;
610          expv_x2 += tks[i].t * tks[i].t * invdt * p;
611          normw += invdt * p;
612        } // setting Z=0 excludes double assignment
613      }
614    }
615
616    mean = mean / normw;
617    expv_x2 = expv_x2 / normw;
618    const double time_var = expv_x2 - mean * mean;
619    const double crappy_error_guess = std::sqrt(time_var);
620    /*GlobalError dummyErrorWithTime(0,
621                                   0,0,
622                                   0,0,0,
623                                   0,0,0,crappy_error_guess);*/
624    //TransientVertex v(pos, time, dummyErrorWithTime, vertexTracks, 5);
625
626    candidate->ClusterIndex = clusterIndex++;
627    ;
628    candidate->Position.SetXYZT(0.0, 0.0, z * 10.0, time * c_light);
629
630    // TBC - fill error later ...
631    candidate->PositionError.SetXYZT(0.0, 0.0, 0.0, crappy_error_guess * c_light);
632
633    clusterIndex++;
634    clusters.push_back(candidate);
635  }
636
637  return clusters;
638}
639
640//------------------------------------------------------------------------------
641
642static double Eik(const track_t &t, const vertex_t &k)
643{
644  return std::pow(t.z - k.z, 2.) / t.dz2 + std::pow(t.t - k.t, 2.) / t.dt2;
645}
646
647//------------------------------------------------------------------------------
648
649static void dump(const double beta, const vector<vertex_t> &y, const vector<track_t> &tks0)
650{
651  // copy and sort for nicer printout
652  vector<track_t> tks;
653  for(vector<track_t>::const_iterator t = tks0.begin(); t != tks0.end(); t++)
654  {
655    tks.push_back(*t);
656  }
657  std::stable_sort(tks.begin(), tks.end(), recTrackLessZ1);
658
659  cout << "-----DAClusterizerInZT::dump ----" << endl;
660  cout << " beta=" << beta << endl;
661  cout << "                                                               z= ";
662  cout.precision(4);
663  for(vector<vertex_t>::const_iterator k = y.begin(); k != y.end(); k++)
664  {
665    //cout  <<  setw(8) << fixed << k->z;
666  }
667  cout << endl
668       << "                                                               t= ";
669  for(vector<vertex_t>::const_iterator k = y.begin(); k != y.end(); k++)
670  {
671    //cout  <<  setw(8) << fixed << k->t;
672  }
673  //cout << endl << "T=" << setw(15) << 1./beta <<"                                             Tc= ";
674  for(vector<vertex_t>::const_iterator k = y.begin(); k != y.end(); k++)
675  {
676    //cout  <<  setw(8) << fixed << k->Tc ;
677  }
678
679  cout << endl
680       << "                                                              pk=";
681  double sumpk = 0;
682  for(vector<vertex_t>::const_iterator k = y.begin(); k != y.end(); k++)
683  {
684    //cout <<  setw(8) <<  setprecision(3) <<  fixed << k->pk;
685    sumpk += k->pk;
686  }
687  cout << endl;
688
689  double E = 0, F = 0;
690  cout << endl;
691  cout << "----       z +/- dz        t +/- dt        ip +/-dip       pt    phi  eta    weights  ----" << endl;
692  cout.precision(4);
693  for(unsigned int i = 0; i < tks.size(); i++)
694  {
695    if(tks[i].Z > 0)
696    {
697      F -= log(tks[i].Z) / beta;
698    }
699    double tz = tks[i].z;
700    double tt = tks[i].t;
701    //cout <<  setw (3)<< i << ")" <<  setw (8) << fixed << setprecision(4)<<  tz << " +/-" <<  setw (6)<< sqrt(tks[i].dz2)
702    //     << setw(8) << fixed << setprecision(4) << tt << " +/-" << setw(6) << std::sqrt(tks[i].dt2)  ;
703
704    double sump = 0.;
705    for(vector<vertex_t>::const_iterator k = y.begin(); k != y.end(); k++)
706    {
707      if((tks[i].pi > 0) && (tks[i].Z > 0))
708      {
709        //double p=pik(beta,tks[i],*k);
710        double p = k->pk * std::exp(-beta * Eik(tks[i], *k)) / tks[i].Z;
711        if(p > 0.0001)
712        {
713          //cout <<  setw (8) <<  setprecision(3) << p;
714        }
715        else
716        {
717          cout << "    .   ";
718        }
719        E += p * Eik(tks[i], *k);
720        sump += p;
721      }
722      else
723      {
724        cout << "        ";
725      }
726    }
727    cout << endl;
728  }
729  cout << endl
730       << "T=" << 1 / beta << " E=" << E << " n=" << y.size() << "  F= " << F << endl
731       << "----------" << endl;
732}
733
734//------------------------------------------------------------------------------
735
736static double update1(double beta, vector<track_t> &tks, vector<vertex_t> &y)
737{
738  //update weights and vertex positions
739  // mass constrained annealing without noise
740  // returns the squared sum of changes of vertex positions
741
742  unsigned int nt = tks.size();
743
744  //initialize sums
745  double sumpi = 0;
746  for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); ++k)
747  {
748    k->sw = 0.;
749    k->swz = 0.;
750    k->swt = 0.;
751    k->se = 0.;
752    k->swE = 0.;
753    k->Tc = 0.;
754  }
755
756  // loop over tracks
757  for(unsigned int i = 0; i < nt; i++)
758  {
759
760    // update pik and Zi
761    double Zi = 0.;
762    for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); ++k)
763    {
764      k->ei = std::exp(-beta * Eik(tks[i], *k)); // cache exponential for one track at a time
765      Zi += k->pk * k->ei;
766    }
767    tks[i].Z = Zi;
768
769    // normalization for pk
770    if(tks[i].Z > 0)
771    {
772      sumpi += tks[i].pi;
773      // accumulate weighted z and weights for vertex update
774      for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); ++k)
775      {
776        k->se += tks[i].pi * k->ei / Zi;
777        const double w = k->pk * tks[i].pi * k->ei / (Zi * (tks[i].dz2 * tks[i].dt2));
778        k->sw += w;
779        k->swz += w * tks[i].z;
780        k->swt += w * tks[i].t;
781        k->swE += w * Eik(tks[i], *k);
782      }
783    }
784    else
785    {
786      sumpi += tks[i].pi;
787    }
788
789  } // end of track loop
790
791  // now update z and pk
792  double delta = 0;
793  for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
794  {
795    if(k->sw > 0)
796    {
797      const double znew = k->swz / k->sw;
798      const double tnew = k->swt / k->sw;
799      delta += std::pow(k->z - znew, 2.) + std::pow(k->t - tnew, 2.);
800      k->z = znew;
801      k->t = tnew;
802      k->Tc = 2. * k->swE / k->sw;
803    }
804    else
805    {
806      // cout << " a cluster melted away ?  pk=" << k->pk <<  " sumw=" << k->sw <<  endl
807      k->Tc = -1;
808    }
809
810    k->pk = k->pk * k->se / sumpi;
811  }
812
813  // return how much the prototypes moved
814  return delta;
815}
816
817//------------------------------------------------------------------------------
818
819static double update2(double beta, vector<track_t> &tks, vector<vertex_t> &y, double &rho0, double dzCutOff)
820{
821  // MVF style, no more vertex weights, update tracks weights and vertex positions, with noise
822  // returns the squared sum of changes of vertex positions
823
824  unsigned int nt = tks.size();
825
826  //initialize sums
827  for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
828  {
829    k->sw = 0.;
830    k->swz = 0.;
831    k->swt = 0.;
832    k->se = 0.;
833    k->swE = 0.;
834    k->Tc = 0.;
835  }
836
837  // loop over tracks
838  for(unsigned int i = 0; i < nt; i++)
839  {
840
841    // update pik and Zi and Ti
842    double Zi = rho0 * std::exp(-beta * (dzCutOff * dzCutOff)); // cut-off (eventually add finite size in time)
843    //double Ti = 0.; // dt0*std::exp(-beta*fDtCutOff);
844    for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
845    {
846      k->ei = std::exp(-beta * Eik(tks[i], *k)); // cache exponential for one track at a time
847      Zi += k->pk * k->ei;
848    }
849    tks[i].Z = Zi;
850
851    // normalization
852    if(tks[i].Z > 0)
853    {
854      // accumulate weighted z and weights for vertex update
855      for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
856      {
857        k->se += tks[i].pi * k->ei / Zi;
858        double w = k->pk * tks[i].pi * k->ei / (Zi * (tks[i].dz2 * tks[i].dt2));
859        k->sw += w;
860        k->swz += w * tks[i].z;
861        k->swt += w * tks[i].t;
862        k->swE += w * Eik(tks[i], *k);
863      }
864    }
865
866  } // end of track loop
867
868  // now update z
869  double delta = 0;
870  for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
871  {
872    if(k->sw > 0)
873    {
874      const double znew = k->swz / k->sw;
875      const double tnew = k->swt / k->sw;
876      delta += std::pow(k->z - znew, 2.) + std::pow(k->t - tnew, 2.);
877      k->z = znew;
878      k->t = tnew;
879      k->Tc = 2 * k->swE / k->sw;
880    }
881    else
882    {
883      // cout << " a cluster melted away ?  pk=" << k->pk <<  " sumw=" << k->sw <<  endl;
884      k->Tc = 0;
885    }
886  }
887
888  // return how much the prototypes moved
889  return delta;
890}
891
892//------------------------------------------------------------------------------
893
894static bool merge(vector<vertex_t> &y)
895{
896  // merge clusters that collapsed or never separated, return true if vertices were merged, false otherwise
897
898  if(y.size() < 2) return false;
899
900  for(vector<vertex_t>::iterator k = y.begin(); (k + 1) != y.end(); k++)
901  {
902    if(std::abs((k + 1)->z - k->z) < 1.e-3 && std::abs((k + 1)->t - k->t) < 1.e-3)
903    { // with fabs if only called after freeze-out (splitAll() at highter T)
904      double rho = k->pk + (k + 1)->pk;
905      if(rho > 0)
906      {
907        k->z = (k->pk * k->z + (k + 1)->z * (k + 1)->pk) / rho;
908        k->t = (k->pk * k->t + (k + 1)->t * (k + 1)->pk) / rho;
909      }
910      else
911      {
912        k->z = 0.5 * (k->z + (k + 1)->z);
913        k->t = 0.5 * (k->t + (k + 1)->t);
914      }
915      k->pk = rho;
916
917      y.erase(k + 1);
918      return true;
919    }
920  }
921
922  return false;
923}
924
925//------------------------------------------------------------------------------
926
927static bool merge(vector<vertex_t> &y, double &beta)
928{
929  // merge clusters that collapsed or never separated,
930  // only merge if the estimated critical temperature of the merged vertex is below the current temperature
931  // return true if vertices were merged, false otherwise
932  if(y.size() < 2) return false;
933
934  for(vector<vertex_t>::iterator k = y.begin(); (k + 1) != y.end(); k++)
935  {
936    if(std::abs((k + 1)->z - k->z) < 2.e-3 && std::abs((k + 1)->t - k->t) < 2.e-3)
937    {
938      double rho = k->pk + (k + 1)->pk;
939      double swE = k->swE + (k + 1)->swE - k->pk * (k + 1)->pk / rho * (std::pow((k + 1)->z - k->z, 2.) + std::pow((k + 1)->t - k->t, 2.));
940      double Tc = 2 * swE / (k->sw + (k + 1)->sw);
941
942      if(Tc * beta < 1)
943      {
944        if(rho > 0)
945        {
946          k->z = (k->pk * k->z + (k + 1)->z * (k + 1)->pk) / rho;
947          k->t = (k->pk * k->t + (k + 1)->t * (k + 1)->pk) / rho;
948        }
949        else
950        {
951          k->z = 0.5 * (k->z + (k + 1)->z);
952          k->t = 0.5 * (k->t + (k + 1)->t);
953        }
954        k->pk = rho;
955        k->sw += (k + 1)->sw;
956        k->swE = swE;
957        k->Tc = Tc;
958        y.erase(k + 1);
959        return true;
960      }
961    }
962  }
963
964  return false;
965}
966
967//------------------------------------------------------------------------------
968
969static bool purge(vector<vertex_t> &y, vector<track_t> &tks, double &rho0, const double beta, const double dzCutOff)
970{
971  // eliminate clusters with only one significant/unique track
972  if(y.size() < 2) return false;
973
974  unsigned int nt = tks.size();
975  double sumpmin = nt;
976  vector<vertex_t>::iterator k0 = y.end();
977  for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
978  {
979    int nUnique = 0;
980    double sump = 0;
981    double pmax = k->pk / (k->pk + rho0 * exp(-beta * dzCutOff * dzCutOff));
982    for(unsigned int i = 0; i < nt; i++)
983    {
984      if(tks[i].Z > 0)
985      {
986        double p = k->pk * std::exp(-beta * Eik(tks[i], *k)) / tks[i].Z;
987        sump += p;
988        if((p > 0.9 * pmax) && (tks[i].pi > 0))
989        {
990          nUnique++;
991        }
992      }
993    }
994
995    if((nUnique < 2) && (sump < sumpmin))
996    {
997      sumpmin = sump;
998      k0 = k;
999    }
1000  }
1001
1002  if(k0 != y.end())
1003  {
1004    //cout << "eliminating prototype at " << k0->z << "," << k0->t << " with sump=" << sumpmin << endl;
1005    //rho0+=k0->pk;
1006    y.erase(k0);
1007    return true;
1008  }
1009  else
1010  {
1011    return false;
1012  }
1013}
1014
1015//------------------------------------------------------------------------------
1016
1017static double beta0(double betamax, vector<track_t> &tks, vector<vertex_t> &y, const double coolingFactor)
1018{
1019
1020  double T0 = 0; // max Tc for beta=0
1021  // estimate critical temperature from beta=0 (T=inf)
1022  unsigned int nt = tks.size();
1023
1024  for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
1025  {
1026
1027    // vertex fit at T=inf
1028    double sumwz = 0.;
1029    double sumwt = 0.;
1030    double sumw = 0.;
1031    for(unsigned int i = 0; i < nt; i++)
1032    {
1033      double w = tks[i].pi / (tks[i].dz2 * tks[i].dt2);
1034      sumwz += w * tks[i].z;
1035      sumwt += w * tks[i].t;
1036      sumw += w;
1037    }
1038    k->z = sumwz / sumw;
1039    k->t = sumwt / sumw;
1040
1041    // estimate Tcrit, eventually do this in the same loop
1042    double a = 0, b = 0;
1043    for(unsigned int i = 0; i < nt; i++)
1044    {
1045      double dx = tks[i].z - (k->z);
1046      double dt = tks[i].t - (k->t);
1047      double w = tks[i].pi / (tks[i].dz2 * tks[i].dt2);
1048      a += w * (std::pow(dx, 2.) / tks[i].dz2 + std::pow(dt, 2.) / tks[i].dt2);
1049      b += w;
1050    }
1051    double Tc = 2. * a / b; // the critical temperature of this vertex
1052    if(Tc > T0) T0 = Tc;
1053  } // vertex loop (normally there should be only one vertex at beta=0)
1054
1055  if(T0 > 1. / betamax)
1056  {
1057    return betamax / pow(coolingFactor, int(std::log(T0 * betamax) / std::log(coolingFactor)) - 1);
1058  }
1059  else
1060  {
1061    // ensure at least one annealing step
1062    return betamax / coolingFactor;
1063  }
1064}
1065
1066//------------------------------------------------------------------------------
1067
1068static bool split(double beta, vector<track_t> &tks, vector<vertex_t> &y)
1069{
1070  // split only critical vertices (Tc >~ T=1/beta   <==>   beta*Tc>~1)
1071  // an update must have been made just before doing this (same beta, no merging)
1072  // returns true if at least one cluster was split
1073
1074  const double epsilon = 1e-3; // split all single vertices by 10 um
1075  bool split = false;
1076
1077  // avoid left-right biases by splitting highest Tc first
1078
1079  std::vector<std::pair<double, unsigned int> > critical;
1080  for(unsigned int ik = 0; ik < y.size(); ik++)
1081  {
1082    if(beta * y[ik].Tc > 1.)
1083    {
1084      critical.push_back(make_pair(y[ik].Tc, ik));
1085    }
1086  }
1087  std::stable_sort(critical.begin(), critical.end(), std::greater<std::pair<double, unsigned int> >());
1088
1089  for(unsigned int ic = 0; ic < critical.size(); ic++)
1090  {
1091    unsigned int ik = critical[ic].second;
1092    // estimate subcluster positions and weight
1093    double p1 = 0, z1 = 0, t1 = 0, w1 = 0;
1094    double p2 = 0, z2 = 0, t2 = 0, w2 = 0;
1095    //double sumpi=0;
1096    for(unsigned int i = 0; i < tks.size(); i++)
1097    {
1098      if(tks[i].Z > 0)
1099      {
1100        //sumpi+=tks[i].pi;
1101        double p = y[ik].pk * exp(-beta * Eik(tks[i], y[ik])) / tks[i].Z * tks[i].pi;
1102        double w = p / (tks[i].dz2 * tks[i].dt2);
1103        if(tks[i].z < y[ik].z)
1104        {
1105          p1 += p;
1106          z1 += w * tks[i].z;
1107          t1 += w * tks[i].t;
1108          w1 += w;
1109        }
1110        else
1111        {
1112          p2 += p;
1113          z2 += w * tks[i].z;
1114          t2 += w * tks[i].t;
1115          w2 += w;
1116        }
1117      }
1118    }
1119    if(w1 > 0)
1120    {
1121      z1 = z1 / w1;
1122      t1 = t1 / w1;
1123    }
1124    else
1125    {
1126      z1 = y[ik].z - epsilon;
1127      t1 = y[ik].t - epsilon;
1128    }
1129    if(w2 > 0)
1130    {
1131      z2 = z2 / w2;
1132      t2 = t2 / w2;
1133    }
1134    else
1135    {
1136      z2 = y[ik].z + epsilon;
1137      t2 = y[ik].t + epsilon;
1138    }
1139
1140    // reduce split size if there is not enough room
1141    if((ik > 0) && (y[ik - 1].z >= z1))
1142    {
1143      z1 = 0.5 * (y[ik].z + y[ik - 1].z);
1144      t1 = 0.5 * (y[ik].t + y[ik - 1].t);
1145    }
1146    if((ik + 1 < y.size()) && (y[ik + 1].z <= z2))
1147    {
1148      z2 = 0.5 * (y[ik].z + y[ik + 1].z);
1149      t2 = 0.5 * (y[ik].t + y[ik + 1].t);
1150    }
1151
1152    // split if the new subclusters are significantly separated
1153    if((z2 - z1) > epsilon || std::abs(t2 - t1) > epsilon)
1154    {
1155      split = true;
1156      vertex_t vnew;
1157      vnew.pk = p1 * y[ik].pk / (p1 + p2);
1158      y[ik].pk = p2 * y[ik].pk / (p1 + p2);
1159      vnew.z = z1;
1160      vnew.t = t1;
1161      y[ik].z = z2;
1162      y[ik].t = t2;
1163      y.insert(y.begin() + ik, vnew);
1164
1165      // adjust remaining pointers
1166      for(unsigned int jc = ic; jc < critical.size(); jc++)
1167      {
1168        if(critical[jc].second > ik)
1169        {
1170          critical[jc].second++;
1171        }
1172      }
1173    }
1174  }
1175
1176  //  stable_sort(y.begin(), y.end(), clusterLessZ);
1177  return split;
1178}
1179
1180//------------------------------------------------------------------------------
1181
1182void splitAll(vector<vertex_t> &y)
1183{
1184
1185  const double epsilon = 1e-3; // split all single vertices by 10 um
1186  const double zsep = 2 * epsilon; // split vertices that are isolated by at least zsep (vertices that haven't collapsed)
1187  const double tsep = 2 * epsilon; // check t as well
1188
1189  vector<vertex_t> y1;
1190
1191  for(vector<vertex_t>::iterator k = y.begin(); k != y.end(); k++)
1192  {
1193    if(((k == y.begin()) || (k - 1)->z < k->z - zsep) && (((k + 1) == y.end()) || (k + 1)->z > k->z + zsep))
1194    {
1195      // isolated prototype, split
1196      vertex_t vnew;
1197      vnew.z = k->z - epsilon;
1198      vnew.t = k->t - epsilon;
1199      (*k).z = k->z + epsilon;
1200      (*k).t = k->t + epsilon;
1201      vnew.pk = 0.5 * (*k).pk;
1202      (*k).pk = 0.5 * (*k).pk;
1203      y1.push_back(vnew);
1204      y1.push_back(*k);
1205    }
1206    else if(y1.empty() || (y1.back().z < k->z - zsep) || (y1.back().t < k->t - tsep))
1207    {
1208      y1.push_back(*k);
1209    }
1210    else
1211    {
1212      y1.back().z -= epsilon;
1213      y1.back().t -= epsilon;
1214      k->z += epsilon;
1215      k->t += epsilon;
1216      y1.push_back(*k);
1217    }
1218  } // vertex loop
1219
1220  y = y1;
1221}
Note: See TracBrowser for help on using the repository browser.