Fork me on GitHub

source: git/modules/VertexFinderDA4D.cc@ 95b4e9f

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 95b4e9f was 95b4e9f, checked in by Pavel Demin <pavel.demin@…>, 8 years ago

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