Fork me on GitHub

source: git/modules/VertexFinder4D.cc@ 3c46e17

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 3c46e17 was 3f8df25, checked in by Michele Selvaggi <michele.selvaggi@…>, 8 years ago

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