Fork me on GitHub

source: git/modules/VertexSorter.cc@ b1fd38d

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

adapt vertexing modules to ROOT 5.34 and GCC 4.4

  • Property mode set to 100644
File size: 7.6 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/ExRootResult.h"
19#include "ExRootAnalysis/ExRootFilter.h"
20#include "ExRootAnalysis/ExRootClassifier.h"
21
22#include "TMath.h"
23#include "TString.h"
24#include "TFormula.h"
25#include "TRandom3.h"
26#include "TObjArray.h"
27#include "TDatabasePDG.h"
28#include "TLorentzVector.h"
29#include "TMatrixT.h"
30#include "TVector3.h"
31
32#include <utility>
33#include <algorithm>
34#include <stdexcept>
35#include <iostream>
36#include <vector>
37#include <map>
38#include <string>
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.