Fork me on GitHub

source: git/modules/VertexFinder.cc@ 5a697dde

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

reorganize includes in vertexing modules
  • Property mode set to 100644
File size: 13.8 KB
Line 
1/** \class VertexFinder
2 *
3 * Cluster vertices from tracks
4 *
5 * \authors A. Hart, M. Selvaggi
6 *
7 */
8
9
10#include "modules/VertexFinder.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#include <map>
36#include <string>
37
38using namespace std;
39
40static const Double_t mm = 1.;
41static const Double_t m = 1000.*mm;
42static const Double_t ns = 1.;
43static const Double_t s = 1.e+9 *ns;
44static const Double_t c_light = 2.99792458e+8 * m/s;
45
46//------------------------------------------------------------------------------
47
48VertexFinder::VertexFinder() :
49 fSigma(0), fMinPT(0), fMaxEta(0), fSeedMinPT(0), fMinNDF(0), fGrowSeeds(0)
50{
51}
52
53//------------------------------------------------------------------------------
54
55VertexFinder::~VertexFinder()
56{
57}
58
59//------------------------------------------------------------------------------
60
61void VertexFinder::Init()
62{
63 fSigma = GetDouble("Sigma", 3.0);
64 fMinPT = GetDouble("MinPT", 0.1);
65 fMaxEta = GetDouble("MaxEta", 10.0);
66 fSeedMinPT = GetDouble("SeedMinPT", 5.0);
67 fMinNDF = GetInt("MinNDF", 4);
68 fGrowSeeds = GetInt("GrowSeeds", 1);
69
70 fInputArray = ImportArray(GetString("InputArray", "TrackSmearing/tracks"));
71 fItInputArray = fInputArray->MakeIterator();
72
73 fOutputArray = ExportArray(GetString("OutputArray", "tracks"));
74 fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
75}
76
77//------------------------------------------------------------------------------
78
79void VertexFinder::Finish()
80{
81 if(fItInputArray) delete fItInputArray;
82}
83
84//------------------------------------------------------------------------------
85
86static Bool_t secondAscending (pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
87{
88 return (pair0.second < pair1.second);
89}
90
91static Bool_t secondDescending (pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
92{
93 return (pair0.second > pair1.second);
94}
95
96//------------------------------------------------------------------------------
97
98void VertexFinder::Process()
99{
100 Candidate *candidate;
101
102 // Clear the track and cluster maps before starting
103 trackIDToDouble.clear ();
104 trackIDToInt.clear ();
105 trackIDToBool.clear ();
106 clusterIDToDouble.clear ();
107 clusterIDToInt.clear ();
108 clusterIDToBool.clear ();
109 trackPT.clear ();
110 clusterSumPT2.clear ();
111
112 // Create the initial cluster seeds
113 createSeeds ();
114
115 // In order of descending seed pt, grow each cluster. If a cluster ends up with
116 // fewer than MinNDF tracks, release the tracks for other clusters to claim.
117 sort (clusterSumPT2.begin (), clusterSumPT2.end (), secondDescending);
118 for (vector<pair<UInt_t, Double_t> >::const_iterator cluster = clusterSumPT2.begin (); cluster != clusterSumPT2.end (); cluster++)
119 {
120 // Skip the cluster if it no longer has any tracks
121 if (!clusterIDToInt.at (cluster->first).at ("ndf"))
122 continue;
123
124 // Grow the cluster if GrowSeeds is true
125 if (fGrowSeeds)
126 growCluster (cluster->first);
127
128 // If the cluster still has fewer than MinNDF tracks, release the tracks;
129 // otherwise, mark the seed track as claimed
130
131 if ((Int_t) clusterIDToInt.at (cluster->first).at ("ndf") < fMinNDF)
132 {
133 for (map<UInt_t, map<string, Int_t> >::iterator track = trackIDToInt.begin (); track != trackIDToInt.end (); track++)
134 {
135 if (track->second.at ("clusterIndex") != (Int_t) cluster->first)
136 continue;
137 track->second["clusterIndex"] = -1;
138 trackIDToBool[track->first]["claimed"] = false;
139 }
140 }
141 else
142 trackIDToBool[clusterIDToInt.at (cluster->first).at ("seed")]["claimed"] = true;
143 }
144
145 // Add tracks to the output array after updating their ClusterIndex.
146 fItInputArray->Reset ();
147 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
148 {
149 if (candidate->Momentum.Pt () < fMinPT || fabs (candidate->Momentum.Eta ()) > fMaxEta)
150 continue;
151 candidate->ClusterIndex = trackIDToInt.at (candidate->GetUniqueID ()).at ("clusterIndex");
152 fOutputArray->Add(candidate);
153 }
154
155 // Add clusters with at least MinNDF tracks to the output array in order of
156 // descending sum(pt**2).
157 clusterSumPT2.clear ();
158 for (map<UInt_t, map<string, Int_t> >::const_iterator cluster = clusterIDToInt.begin (); cluster != clusterIDToInt.end (); cluster++)
159 {
160
161 if (cluster->second.at ("ndf") < fMinNDF)
162 continue;
163 clusterSumPT2.push_back (make_pair (cluster->first, clusterIDToDouble.at (cluster->first).at ("sumPT2")));
164 }
165 sort (clusterSumPT2.begin (), clusterSumPT2.end (), secondDescending);
166
167 for (vector<pair<UInt_t, Double_t> >::const_iterator cluster = clusterSumPT2.begin (); cluster != clusterSumPT2.end (); cluster++)
168 {
169 DelphesFactory *factory = GetFactory();
170 candidate = factory->NewCandidate();
171
172 candidate->ClusterIndex = cluster->first;
173 candidate->ClusterNDF = clusterIDToInt.at (cluster->first).at ("ndf");
174 candidate->ClusterSigma = fSigma;
175 candidate->SumPT2 = cluster->second;
176 candidate->Position.SetXYZT(0.0, 0.0, clusterIDToDouble.at (cluster->first).at ("z"), 0.0);
177 candidate->PositionError.SetXYZT(0.0, 0.0, clusterIDToDouble.at (cluster->first).at ("ez"), 0.0);
178
179 fVertexOutputArray->Add(candidate);
180 }
181}
182
183//------------------------------------------------------------------------------
184
185void VertexFinder::createSeeds ()
186{
187 Candidate *candidate;
188 UInt_t clusterIndex = 0, maxSeeds = 0;
189
190 // Loop over all tracks, initializing some variables.
191 fItInputArray->Reset();
192 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
193 {
194 if (candidate->Momentum.Pt () < fMinPT || fabs (candidate->Momentum.Eta ()) > fMaxEta)
195 continue;
196
197 trackIDToDouble[candidate->GetUniqueID ()]["pt"] = candidate->Momentum.Pt ();
198 trackIDToDouble[candidate->GetUniqueID ()]["ept"] = candidate->ErrorPT ? candidate->ErrorPT : 1.0e-15;;
199 trackIDToDouble[candidate->GetUniqueID ()]["eta"] = candidate->Momentum.Eta ();
200
201 trackIDToDouble[candidate->GetUniqueID ()]["z"] = candidate->DZ;
202 trackIDToDouble[candidate->GetUniqueID ()]["ez"] = candidate->ErrorDZ ? candidate->ErrorDZ : 1.0e-15;
203
204 trackIDToInt[candidate->GetUniqueID ()]["clusterIndex"] = -1;
205 trackIDToInt[candidate->GetUniqueID ()]["interactionIndex"] = candidate->IsPU;
206
207 trackIDToBool[candidate->GetUniqueID ()]["claimed"] = false;
208
209 trackPT.push_back (make_pair (candidate->GetUniqueID (), candidate->Momentum.Pt ()));
210 }
211
212 // Sort tracks by pt and leave only the SeedMinPT highest pt ones in the
213 // trackPT vector.
214 sort (trackPT.begin (), trackPT.end (), secondDescending);
215 for (vector<pair<UInt_t, Double_t> >::const_iterator track = trackPT.begin (); track != trackPT.end (); track++, maxSeeds++)
216 {
217 if (track->second < fSeedMinPT)
218 break;
219 }
220 // If there are no tracks with pt above MinSeedPT, create just one seed from
221 // the highest pt track.
222 if (!maxSeeds)
223 maxSeeds++;
224 if (trackPT.size () > maxSeeds)
225 {
226 trackPT.erase (trackPT.begin () + maxSeeds, trackPT.end ());
227 }
228
229 // Create the seeds from the SeedMinPT highest pt tracks.
230 for (vector<pair<UInt_t, Double_t> >::const_iterator track = trackPT.begin (); track != trackPT.end (); track++, clusterIndex++)
231 {
232 addTrackToCluster (track->first, clusterIndex);
233 clusterSumPT2.push_back (make_pair (clusterIndex, track->second * track->second));
234 }
235}
236
237//------------------------------------------------------------------------------
238
239void VertexFinder::growCluster (const UInt_t clusterIndex)
240{
241 Bool_t done = false;
242 UInt_t nearestID;
243 Int_t oldClusterIndex;
244 Double_t nearestDistance;
245 vector<UInt_t> nearTracks;
246 nearTracks.clear ();
247
248 // Grow the cluster until there are no more tracks within Sigma standard
249 // deviations of the cluster.
250 while (!done)
251 {
252 done = true;
253 nearestID = 0;
254 nearestDistance = -1.0;
255
256 // These two loops are for finding the nearest track to the cluster. The
257 // first time, the ID of each track within 10*Sigma of the cluster is
258 // saved in the nearTracks vector; subsequently, to save time, only the
259 // tracks in this vector are checked.
260 if (!nearTracks.size ())
261 {
262
263 for (map<UInt_t, map<string, Double_t> >::const_iterator track = trackIDToDouble.begin (); track != trackIDToDouble.end (); track++)
264 {
265 if (trackIDToBool.at (track->first).at ("claimed") || trackIDToInt.at (track->first).at ("clusterIndex") == (Int_t) clusterIndex)
266 continue;
267
268 Double_t distance = fabs (clusterIDToDouble.at (clusterIndex).at ("z") - track->second.at ("z")) / hypot (clusterIDToDouble.at (clusterIndex).at ("ez"), track->second.at ("ez"));
269 if (nearestDistance < 0.0 || distance < nearestDistance)
270 {
271 nearestID = track->first;
272 nearestDistance = distance;
273 }
274 if (distance < 10.0 * fSigma)
275 nearTracks.push_back (track->first);
276 }
277 }
278
279 else
280 {
281 for (vector<UInt_t>::const_iterator track = nearTracks.begin (); track != nearTracks.end (); track++)
282 {
283 if (trackIDToBool.at (*track).at ("claimed") || trackIDToInt.at (*track).at ("clusterIndex") == (Int_t) clusterIndex)
284 continue;
285 Double_t distance = fabs (clusterIDToDouble.at (clusterIndex).at ("z") - trackIDToDouble.at (*track).at ("z")) / hypot (clusterIDToDouble.at (clusterIndex).at ("ez"), trackIDToDouble.at (*track).at ("ez"));
286 if (nearestDistance < 0.0 || distance < nearestDistance)
287 {
288 nearestID = *track;
289 nearestDistance = distance;
290 }
291 }
292 }
293
294 // If no tracks within Sigma of the cluster were found, stop growing.
295 done = nearestDistance > fSigma || nearestDistance < 0.0;
296 if (done)
297 {
298 continue;
299 }
300
301 // Add the nearest track within Sigma to the cluster. If it already
302 // belonged to another cluster, remove it from that cluster first.
303 if (nearestDistance < fSigma)
304 {
305 oldClusterIndex = trackIDToInt.at (nearestID).at ("clusterIndex");
306 if (oldClusterIndex >= 0)
307 removeTrackFromCluster (nearestID, oldClusterIndex);
308
309 trackIDToBool[nearestID]["claimed"] = true;
310 addTrackToCluster (nearestID, clusterIndex);
311 }
312 }
313}
314
315//------------------------------------------------------------------------------
316
317Double_t VertexFinder::weight (const UInt_t trackID)
318{
319 return ((trackIDToDouble.at (trackID).at ("pt") / (trackIDToDouble.at (trackID).at ("ept") * trackIDToDouble.at (trackID).at ("ez"))) * (trackIDToDouble.at (trackID).at ("pt") / (trackIDToDouble.at (trackID).at ("ept") * trackIDToDouble.at (trackID).at ("ez"))));
320}
321
322//------------------------------------------------------------------------------
323
324void VertexFinder::removeTrackFromCluster (const UInt_t trackID, const UInt_t clusterID)
325{
326 Double_t wz = weight (trackID);
327
328 trackIDToInt[trackID]["clusterIndex"] = -1;
329 clusterIDToInt[clusterID]["ndf"]--;
330
331 clusterIDToDouble[clusterID]["sumZ"] -= wz * trackIDToDouble.at (trackID).at ("z");
332 clusterIDToDouble[clusterID]["errorSumZ"] -= wz * trackIDToDouble.at (trackID).at ("ez") * trackIDToDouble.at (trackID).at ("ez");
333 clusterIDToDouble[clusterID]["sumOfWeightsZ"] -= wz;
334 clusterIDToDouble[clusterID]["z"] = clusterIDToDouble.at (clusterID).at ("sumZ") / clusterIDToDouble.at (clusterID).at ("sumOfWeightsZ");
335 clusterIDToDouble[clusterID]["ez"] = sqrt ((1.0 / clusterIDToInt.at (clusterID).at ("ndf")) * (clusterIDToDouble.at (clusterID).at ("errorSumZ") / clusterIDToDouble.at (clusterID).at ("sumOfWeightsZ")));
336 clusterIDToDouble[clusterID]["sumPT2"] -= trackIDToDouble.at (trackID).at ("pt") * trackIDToDouble.at (trackID).at ("pt");
337}
338
339//------------------------------------------------------------------------------
340
341void VertexFinder::addTrackToCluster (const UInt_t trackID, const UInt_t clusterID)
342{
343 Double_t wz = weight (trackID);
344
345 if (!clusterIDToInt.count (clusterID))
346 {
347 clusterIDToInt[clusterID]["ndf"] = 0;
348 clusterIDToInt[clusterID]["seed"] = trackID;
349 clusterIDToDouble[clusterID]["sumZ"] = 0.0;
350 clusterIDToDouble[clusterID]["errorSumZ"] = 0.0;
351 clusterIDToDouble[clusterID]["sumOfWeightsZ"] = 0.0;
352 clusterIDToDouble[clusterID]["sumPT2"] = 0.0;
353 }
354
355 trackIDToInt[trackID]["clusterIndex"] = clusterID;
356 clusterIDToInt[clusterID]["ndf"]++;
357
358 clusterIDToDouble[clusterID]["sumZ"] += wz * trackIDToDouble.at (trackID).at ("z");
359 clusterIDToDouble[clusterID]["errorSumZ"] += wz * trackIDToDouble.at (trackID).at ("ez") * trackIDToDouble.at (trackID).at ("ez");
360 clusterIDToDouble[clusterID]["sumOfWeightsZ"] += wz;
361 clusterIDToDouble[clusterID]["z"] = clusterIDToDouble.at (clusterID).at ("sumZ") / clusterIDToDouble.at (clusterID).at ("sumOfWeightsZ");
362 clusterIDToDouble[clusterID]["ez"] = sqrt ((1.0 / clusterIDToInt.at (clusterID).at ("ndf")) * (clusterIDToDouble.at (clusterID).at ("errorSumZ") / clusterIDToDouble.at (clusterID).at ("sumOfWeightsZ")));
363 clusterIDToDouble[clusterID]["sumPT2"] += trackIDToDouble.at (trackID).at ("pt") * trackIDToDouble.at (trackID).at ("pt");
364}
365
366//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.