Fork me on GitHub

source: git/modules/VertexFinderDA4D.h@ 77249aa

Timing
Last change on this file since 77249aa was 4ac0049, checked in by Kaan Yüksel Oyulmaz <kaanyukseloyulmaz@…>, 5 years ago

VertexFinderDA4D.cc module is cleaned from graphs, FCChh_PileUpVtx.tcl card, PileUpSubtractor4D.cc and TimeSmearing.cc were updated in meeting

  • Property mode set to 100644
File size: 7.5 KB
Line 
1#ifndef VertexFinderDA4D_h
2#define VertexFinderDA4D_h
3
4/** \class VertexFinderDA4D
5 *
6 * Cluster vertices from tracks using deterministic annealing and timing information
7 *
8 * Author O. Cerri
9 *
10 */
11
12#include "classes/DelphesModule.h"
13
14#include <vector>
15#include <iostream>
16
17class TObjArray;
18class TIterator;
19class Candidate;
20
21class VertexFinderDA4D: public DelphesModule
22{
23 public:
24
25 class tracks_t
26 {
27 public:
28 std::vector<double> z; // z-coordinate at point of closest approach to the beamline [mm]
29 std::vector<double> r; // z-coordinate at point of closest approach to the beamline [mm]
30
31 std::vector<double> t; // t-coordinate at point of closest approach to the beamline [ps]
32
33 std::vector<double> dz2_o; //1 over the squared error of z(pca)
34 std::vector<double> dz_o;
35 std::vector<double> dt2_o; //1 over the squared error of t(pca)
36 std::vector<double> dt_o;
37
38 std::vector<Candidate*> tt; // a pointer to the Candidate Track
39 std::vector<double> Z; // Normalization of P(y|x) = p_k * exp(-beta*E_ik)/Z_i
40
41 std::vector<double> w; // track weight
42
43 // std::vector<double> pt;
44 // std::vector<double> pz;
45
46 std::vector<int> PID;
47
48 double sum_w_o_dz2 = 0;
49 double sum_w_o_dt2 = 0;
50 double sum_w = 0;
51
52 tracks_t(){}
53 ~tracks_t(){}
54
55 void addItem( double new_z, double new_t, double new_dz2_o, double new_dt2_o, Candidate * new_tt, double new_w, int new_PID)
56 {
57 z.push_back( new_z );
58 t.push_back( new_t );
59 dz2_o.push_back( new_dz2_o );
60 dz_o.push_back( sqrt(new_dz2_o) );
61 dt2_o.push_back( new_dt2_o );
62 dt_o.push_back( sqrt(new_dt2_o) );
63 tt.push_back( new_tt );
64
65 w.push_back( new_w );
66 Z.push_back( 1.0 );
67
68 PID.push_back(new_PID);
69 }
70
71 unsigned int getSize()
72 {
73 return z.size();
74 }
75
76 };
77
78 class vertex_t
79 {
80 public:
81 std::vector<double> z; // z coordinate
82 std::vector<double> t; // t coordinate
83 std::vector<double> pk; // vertex weight for "constrained" clustering
84
85 // Elements of the covariance matrix of the posterior distribution
86 std::vector<double> szz, stt, stz;
87 std::vector<double> beta_c;
88
89 double ZSize;
90 double TSize;
91
92 // Used in the tracks-cluster assignment
93 double sum_pt;
94 std::vector<unsigned int> i_tracks;
95
96
97 vertex_t(){}
98 ~vertex_t(){}
99
100 void addItem( double new_z, double new_t, double new_pk)
101 {
102 z.push_back( new_z);
103 t.push_back( new_t);
104 pk.push_back( new_pk);
105
106 szz.push_back(0);
107 stt.push_back(0);
108 stz.push_back(0);
109 beta_c.push_back(0);
110 }
111
112 void removeItem( unsigned int i )
113 {
114 z.erase( z.begin() + i );
115 t.erase( t.begin() + i );
116 pk.erase( pk.begin() + i );
117
118 szz.erase(szz.begin() + i);
119 stt.erase(stt.begin() + i);
120 stz.erase(stz.begin() + i);
121 beta_c.erase(beta_c.begin() + i);
122 }
123
124 void mergeItems(unsigned int k1, unsigned int k2)
125 {
126 double new_pk = pk[k1] + pk[k2];
127 double new_z = (z[k1]*pk[k1] + z[k2]*pk[k2])/new_pk;
128 double new_t = (t[k1]*pk[k1] + t[k2]*pk[k2])/new_pk;
129
130 removeItem(k2);
131 z[k1] = new_z;
132 t[k1] = new_t;
133 pk[k1] = new_pk;
134 }
135
136 unsigned int getSize() const
137 {
138 return z.size();
139 }
140
141 double DistanceSquare(unsigned int k1, unsigned int k2)
142 {
143 double dz = (z[k1] - z[k2])/ZSize;
144 double dt = (t[k1] - t[k2])/TSize;
145
146 return dz*dz + dt*dt;
147 }
148
149 unsigned int NearestCluster(double t_, double z_)
150 {
151 unsigned int k_min = 0;
152 double d2_min = 0;
153
154 for(unsigned int k = 0; k < getSize(); k++)
155 {
156 double dt = (t[k] - t_)/TSize;
157 double dz = (z[k] - z_)/ZSize;
158 double d2 = dt*dt + dz*dz;
159 if(k == 0 || d2 < d2_min)
160 {
161 k_min = k;
162 d2_min = d2;
163 }
164 }
165
166 return k_min;
167 }
168
169 std::pair<double, unsigned int> ComputeAllBeta_c(unsigned int verbose = 0)
170 {
171 unsigned int nv = getSize();
172 unsigned int k_min = 0;
173 double beta_c_min = 0;
174 for(unsigned int k = 0; k < nv; k++)
175 {
176 if(szz[k] == 0 || stt[k] == 0)
177 {
178 beta_c[k] = 1000;
179 if(verbose>5)
180 {
181 std::cout << std::endl << Form("Varianca too small! beta_c set to 1000\nVertex %d: t=%1.2e, z=%1.2e, pk=%1.2e, szz=%1.2e, stt=%1.2e, stz=%1.2e", k, t[k], z[k], pk[k], szz[k], stt[k], stz[k]);
182 // throw std::invalid_argument( "Attempting to compute beta c for uninitialized vertex" );
183 }
184 }
185
186 double aux = (stt[k] - szz[k])*(stt[k] - szz[k]) + 4*stz[k]*stz[k];
187 aux = 1. / (stt[k] + szz[k] + sqrt(aux));
188
189 if(k == 0 || aux < beta_c_min)
190 {
191 beta_c_min = aux;
192 k_min = k;
193 }
194
195 beta_c[k] = aux;
196 }
197
198 std::pair<double, unsigned int> out(beta_c_min, k_min);
199 return out;
200 }
201 };
202
203 VertexFinderDA4D();
204 ~VertexFinderDA4D();
205
206 void Init();
207 void Process();
208 void Finish();
209
210 void clusterize(TObjArray &clusters);
211
212 // Define the distance metric between tracks and vertices
213 double Energy(double t_z, double v_z, double dz2_o, double t_t, double v_t, double dt2_o);
214
215 // Fill the tracks structure from the input array
216 void fill(tracks_t &tks);
217
218 // Compute higher phase transition temperature
219 double beta0(tracks_t & tks, vertex_t &vtx);
220
221 // Compute the new vertexes position and mass
222 double update(double beta, tracks_t &tks, vertex_t &vtx, double rho0);
223
224 // If a vertex has beta_c lower than beta, split it
225 bool split(double &beta, vertex_t & vtx, tracks_t & tks);
226
227 // Merge vertexes closer than declared dimensions
228 bool merge(vertex_t & vtx, double d2_merge);
229
230 // Eliminate clusters with only one significant/unique track
231 bool purge(vertex_t & vtx, tracks_t & tks, double & rho0, const double beta, double min_prob, double min_trk);
232
233 // Compute all the energies and set the partition function normalization for each track
234 std::vector<double> Compute_pk_exp_mBetaE(double beta, vertex_t &vtx, tracks_t &tks, double Z_init);
235
236
237 private:
238
239 unsigned int fVerbose;
240
241 UInt_t fMaxIterations;
242 UInt_t fMaxVertexNumber;
243
244 Float_t fBetaMax;
245 Float_t fBetaStop;
246 Float_t fBetaPurge;
247
248
249 Float_t fVertexZSize;
250 Float_t fVertexTSize;
251
252 Double_t fCoolingFactor;
253
254 Double_t fDzCutOff;
255 Double_t fD0CutOff;
256 Double_t fDtCutOff;
257 Double_t fPtMin;
258 Double_t fPtMax;
259
260 Double_t fD2UpdateLim;
261 Double_t fD2Merge;
262 Double_t fMuOutlayer;
263 Double_t fMinTrackProb;
264 Int_t fMinNTrack;
265
266 TObjArray *fInputArray;
267 TIterator *fItInputArray;
268
269 TObjArray *fTrackOutputArray;
270 TObjArray *fVertexOutputArray;
271
272 ClassDef(VertexFinderDA4D, 1)
273
274 //Used to keep track of number of verices in verbose mode
275 std::vector<double> fEnergy_rec;
276 std::vector<double> fBeta_rec;
277 std::vector<double> fNvtx_rec;
278
279 // Debug purpose only
280 TObjArray *fInputGenVtx;
281 TIterator *fItInputGenVtx;
282
283 TString fFigFolderPath = "";
284};
285
286#endif
Note: See TracBrowser for help on using the repository browser.