Fork me on GitHub

source: git/modules/VertexFinderDA4D.cc@ 3cfe61d

Last change on this file since 3cfe61d was 77e9ae1, checked in by Pavel Demin <pavel-demin@…>, 6 years 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.