Fork me on GitHub

source: git/modules/VertexSorter.cc@ 414db83

ImprovedOutputFile Timing
Last change on this file since 414db83 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: 7.2 KB
Line 
1/** \class VertexSorter
2 *
3 *
4 * Sorts vertices according to different criteria
5 *
6 * \authors A. Hart, M. Selvaggi
7 *
8 *
9*/
10
11#include "modules/VertexSorter.h"
12
13#include "classes/DelphesClasses.h"
14#include "classes/DelphesFactory.h"
15#include "classes/DelphesFormula.h"
16#include "classes/DelphesPileUpReader.h"
17
18#include "ExRootAnalysis/ExRootClassifier.h"
19#include "ExRootAnalysis/ExRootFilter.h"
20#include "ExRootAnalysis/ExRootResult.h"
21
22#include "TDatabasePDG.h"
23#include "TFormula.h"
24#include "TLorentzVector.h"
25#include "TMath.h"
26#include "TMatrixT.h"
27#include "TObjArray.h"
28#include "TRandom3.h"
29#include "TString.h"
30#include "TVector3.h"
31
32#include <algorithm>
33#include <iostream>
34#include <map>
35#include <stdexcept>
36#include <string>
37#include <utility>
38#include <vector>
39
40using namespace std;
41
42static const Double_t mm = 1.;
43static const Double_t m = 1000. * mm;
44static const Double_t ns = 1.;
45static const Double_t s = 1.e+9 * ns;
46static const Double_t c_light = 2.99792458e+8 * m / s;
47
48//------------------------------------------------------------------------------
49
50VertexSorter::VertexSorter() :
51 fInputArray(NULL), fTrackInputArray(NULL), fItTrackInputArray(NULL), fJetInputArray(NULL), fItJetInputArray(NULL), fOutputArray(NULL)
52{
53}
54
55//------------------------------------------------------------------------------
56
57VertexSorter::~VertexSorter()
58{
59}
60
61//------------------------------------------------------------------------------
62
63void VertexSorter::Init()
64{
65 fInputArray = ImportArray(GetString("InputArray", "VertexFinder/vertices"));
66
67 fTrackInputArray = ImportArray(GetString("TrackInputArray", "VertexFinder/tracks"));
68 fItTrackInputArray = fTrackInputArray->MakeIterator();
69
70 if(string(GetString("JetInputArray", "")) != "")
71 {
72 fJetInputArray = ImportArray(GetString("JetInputArray", ""));
73 fItJetInputArray = fJetInputArray->MakeIterator();
74 }
75
76 // import beamspot
77 try
78 {
79 fBeamSpotInputArray = ImportArray(GetString("BeamSpotInputArray", "BeamSpotFilter/beamSpotParticle"));
80 }
81 catch(runtime_error &e)
82 {
83 fBeamSpotInputArray = 0;
84 }
85
86 fOutputArray = ExportArray(GetString("OutputArray", "clusters"));
87
88 fMethod = GetString("Method", "BTV");
89}
90
91//------------------------------------------------------------------------------
92
93void VertexSorter::Finish()
94{
95 if(fItTrackInputArray) delete fItTrackInputArray;
96 if(fItJetInputArray) delete fItJetInputArray;
97}
98
99//------------------------------------------------------------------------------
100
101static Bool_t secondDescending(pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
102{
103 return (pair0.second > pair1.second);
104}
105
106static Bool_t secondAscending(pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
107{
108 return (pair0.second < pair1.second);
109}
110
111//------------------------------------------------------------------------------
112
113void VertexSorter::Process()
114{
115 Candidate *candidate, *jetCandidate, *beamSpotCandidate;
116 map<Int_t, UInt_t> clusterIDToIndex;
117 map<Int_t, UInt_t>::const_iterator itClusterIDToIndex;
118 map<Int_t, Double_t> clusterIDToSumPT2;
119 map<Int_t, Double_t>::const_iterator itClusterIDToSumPT2;
120 vector<pair<Int_t, Double_t> > sortedClusterIDs;
121 vector<pair<Int_t, Double_t> >::const_iterator itSortedClusterIDs;
122
123 for(Int_t iCluster = 0; iCluster < fInputArray->GetEntries(); iCluster++)
124 {
125 const Candidate &cluster = *((Candidate *)fInputArray->At(iCluster));
126 clusterIDToIndex[cluster.ClusterIndex] = iCluster;
127 clusterIDToSumPT2[cluster.ClusterIndex] = 0.0;
128 }
129
130 if(fMethod == "BTV")
131 {
132 if(!fJetInputArray)
133 {
134 cout << "BTV PV sorting selected, but no jet collection given!" << endl;
135 throw 0;
136 }
137
138 fItTrackInputArray->Reset();
139 while((candidate = static_cast<Candidate *>(fItTrackInputArray->Next())))
140 {
141 if(candidate->Momentum.Pt() < 1.0)
142 continue;
143 if(candidate->ClusterIndex < 0)
144 continue;
145 TLorentzVector p(candidate->Momentum.Px(), candidate->Momentum.Py(), candidate->Momentum.Pz(), candidate->Momentum.E());
146 Bool_t isInJet = false;
147
148 fItJetInputArray->Reset();
149 while((jetCandidate = static_cast<Candidate *>(fItJetInputArray->Next())))
150 {
151 if(jetCandidate->Momentum.Pt() < 30.0)
152 continue;
153 TLorentzVector q(jetCandidate->Momentum.Px(), jetCandidate->Momentum.Py(), jetCandidate->Momentum.Pz(), jetCandidate->Momentum.E());
154
155 if(p.DeltaR(q) > 0.4)
156 continue;
157 isInJet = true;
158 break;
159 }
160 if(!isInJet)
161 continue;
162
163 clusterIDToSumPT2.at(candidate->ClusterIndex) += candidate->Momentum.Pt() * candidate->Momentum.Pt();
164 }
165
166 for(itClusterIDToSumPT2 = clusterIDToSumPT2.begin(); itClusterIDToSumPT2 != clusterIDToSumPT2.end(); ++itClusterIDToSumPT2)
167 sortedClusterIDs.push_back(make_pair(itClusterIDToSumPT2->first, itClusterIDToSumPT2->second));
168 sort(sortedClusterIDs.begin(), sortedClusterIDs.end(), secondDescending);
169 }
170 else if(fMethod == "GenClosest")
171 {
172 if(!fBeamSpotInputArray)
173 {
174 cout << "GenClosest PV sorting selected, but no beamspot collection given!" << endl;
175 throw 0;
176 }
177 if(!fBeamSpotInputArray->GetSize())
178 {
179 cout << "Beamspot collection is empty!" << endl;
180 throw 0;
181 }
182
183 beamSpotCandidate = (Candidate *)fBeamSpotInputArray->At(0);
184 for(Int_t iCluster = 0; iCluster < fInputArray->GetEntries(); iCluster++)
185 {
186 const Candidate &cluster = *((Candidate *)fInputArray->At(iCluster));
187 sortedClusterIDs.push_back(make_pair(cluster.ClusterIndex, fabs(cluster.Position.Z() - beamSpotCandidate->Position.Z())));
188 }
189 sort(sortedClusterIDs.begin(), sortedClusterIDs.end(), secondAscending);
190 }
191 else if(fMethod == "GenBest")
192 {
193 fItTrackInputArray->Reset();
194 while((candidate = static_cast<Candidate *>(fItTrackInputArray->Next())))
195 {
196 if(candidate->IsPU)
197 continue;
198 for(itClusterIDToIndex = clusterIDToIndex.begin(); itClusterIDToIndex != clusterIDToIndex.end(); ++itClusterIDToIndex)
199 {
200 if(candidate->ClusterIndex != itClusterIDToIndex->first)
201 continue;
202 clusterIDToSumPT2.at(itClusterIDToIndex->first) += candidate->Momentum.Pt() * candidate->Momentum.Pt();
203 }
204 }
205
206 for(itClusterIDToSumPT2 = clusterIDToSumPT2.begin(); itClusterIDToSumPT2 != clusterIDToSumPT2.end(); ++itClusterIDToSumPT2)
207 sortedClusterIDs.push_back(make_pair(itClusterIDToSumPT2->first, itClusterIDToSumPT2->second));
208 sort(sortedClusterIDs.begin(), sortedClusterIDs.end(), secondDescending);
209 }
210 else
211 {
212 cout << "\"" << fMethod << "\" is not a valid sorting method!" << endl;
213 cout << "Valid methods are:" << endl;
214 cout << " BTV" << endl;
215 cout << " GenClosest" << endl;
216 cout << " GenBest" << endl;
217 throw 0;
218 }
219 for(itSortedClusterIDs = sortedClusterIDs.begin(); itSortedClusterIDs != sortedClusterIDs.end(); ++itSortedClusterIDs)
220 {
221 Candidate *cluster = (Candidate *)fInputArray->At(clusterIDToIndex.at(itSortedClusterIDs->first));
222 if(fMethod == "BTV")
223 cluster->BTVSumPT2 = itSortedClusterIDs->second;
224 else if(fMethod == "GenClosest")
225 cluster->GenDeltaZ = itSortedClusterIDs->second;
226 else if(fMethod == "GenBest")
227 cluster->GenSumPT2 = itSortedClusterIDs->second;
228 fOutputArray->Add(cluster);
229 }
230}
231
232//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.