Fork me on GitHub

source: git/modules/VertexSorter.cc@ 8f0b34c

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 8f0b34c was 2118a6a, checked in by Michele Selvaggi <michele.selvaggi@…>, 9 years ago

fixed missing vertexing variables
  • Property mode set to 100644
File size: 6.9 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
32static const Double_t mm = 1.;
33static const Double_t m = 1000.*mm;
34static const Double_t ns = 1.;
35static const Double_t s = 1.e+9 *ns;
36static const Double_t c_light = 2.99792458e+8 * m/s;
37
38//------------------------------------------------------------------------------
39
40VertexSorter::VertexSorter() :
41 fInputArray(NULL), fTrackInputArray(NULL), fItTrackInputArray(NULL), fJetInputArray(NULL), fItJetInputArray(NULL), fOutputArray(NULL)
42{
43}
44
45//------------------------------------------------------------------------------
46
47VertexSorter::~VertexSorter()
48{
49}
50
51//------------------------------------------------------------------------------
52
53void VertexSorter::Init()
54{
55 fInputArray = ImportArray(GetString("InputArray", "VertexFinder/vertices"));
56
57 fTrackInputArray = ImportArray(GetString("TrackInputArray", "VertexFinder/tracks"));
58 fItTrackInputArray = fTrackInputArray->MakeIterator();
59
60 if (string (GetString("JetInputArray", "")) != "")
61 {
62 fJetInputArray = ImportArray(GetString("JetInputArray", ""));
63 fItJetInputArray = fJetInputArray->MakeIterator();
64 }
65
66 // import beamspot
67 try
68 {
69 fBeamSpotInputArray = ImportArray(GetString("BeamSpotInputArray", "BeamSpotFilter/beamSpotParticle"));
70 }
71 catch(runtime_error &e)
72 {
73 fBeamSpotInputArray = 0;
74 }
75
76 fOutputArray = ExportArray(GetString("OutputArray", "clusters"));
77
78 fMethod = GetString ("Method", "BTV");
79}
80
81//------------------------------------------------------------------------------
82
83void VertexSorter::Finish()
84{
85 if(fItTrackInputArray) delete fItTrackInputArray;
86 if(fItJetInputArray) delete fItJetInputArray;
87}
88
89//------------------------------------------------------------------------------
90//
91Bool_t VertexSorter::secondDescending (pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
92{
93 return (pair0.second > pair1.second);
94}
95Bool_t VertexSorter::secondAscending (pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
96{
97 return (pair0.second < pair1.second);
98}
99
100void VertexSorter::Process()
101{
102 Candidate *candidate, *jetCandidate, *beamSpotCandidate;
103 unordered_map<Int_t, UInt_t> clusterIDToIndex;
104 unordered_map<Int_t, Double_t> clusterIDToSumPT2;
105 vector<pair<Int_t, Double_t> > sortedClusterIDs;
106
107 for (Int_t iCluster = 0; iCluster < fInputArray->GetEntries (); iCluster++)
108 {
109 const Candidate &cluster = *((Candidate *) fInputArray->At (iCluster));
110 clusterIDToIndex[cluster.ClusterIndex] = iCluster;
111 clusterIDToSumPT2[cluster.ClusterIndex] = 0.0;
112 }
113
114 if (fMethod == "BTV")
115 {
116 if (!fJetInputArray)
117 {
118 cout << "BTV PV sorting selected, but no jet collection given!" << endl;
119 throw 0;
120 }
121
122 fItTrackInputArray->Reset();
123 while((candidate = static_cast<Candidate*>(fItTrackInputArray->Next())))
124 {
125 if (candidate->Momentum.Pt () < 1.0)
126 continue;
127 if (candidate->ClusterIndex < 0)
128 continue;
129 TLorentzVector p (candidate->Momentum.Px (), candidate->Momentum.Py (), candidate->Momentum.Pz (), candidate->Momentum.E ());
130 Bool_t isInJet = false;
131
132 fItJetInputArray->Reset();
133 while((jetCandidate = static_cast<Candidate*>(fItJetInputArray->Next())))
134 {
135 if (jetCandidate->Momentum.Pt () < 30.0)
136 continue;
137 TLorentzVector q (jetCandidate->Momentum.Px (), jetCandidate->Momentum.Py (), jetCandidate->Momentum.Pz (), jetCandidate->Momentum.E ());
138
139 if (p.DeltaR (q) > 0.4)
140 continue;
141 isInJet = true;
142 break;
143 }
144 if (!isInJet)
145 continue;
146
147 clusterIDToSumPT2.at (candidate->ClusterIndex) += candidate->Momentum.Pt () * candidate->Momentum.Pt ();
148 }
149
150 for (const auto &clusterID : clusterIDToSumPT2)
151 sortedClusterIDs.push_back (make_pair (clusterID.first, clusterID.second));
152 sort (sortedClusterIDs.begin (), sortedClusterIDs.end (), secondDescending);
153 }
154 else if (fMethod == "GenClosest")
155 {
156 if (!fBeamSpotInputArray)
157 {
158 cout << "GenClosest PV sorting selected, but no beamspot collection given!" << endl;
159 throw 0;
160 }
161 if (!fBeamSpotInputArray->GetSize ())
162 {
163 cout << "Beamspot collection is empty!" << endl;
164 throw 0;
165 }
166
167 beamSpotCandidate = (Candidate *) fBeamSpotInputArray->At (0);
168 for (Int_t iCluster = 0; iCluster < fInputArray->GetEntries (); iCluster++)
169 {
170 const Candidate &cluster = *((Candidate *) fInputArray->At (iCluster));
171 sortedClusterIDs.push_back (make_pair (cluster.ClusterIndex, fabs (cluster.Position.Z () - beamSpotCandidate->Position.Z ())));
172 }
173 sort (sortedClusterIDs.begin (), sortedClusterIDs.end (), secondAscending);
174 }
175 else if (fMethod == "GenBest")
176 {
177 fItTrackInputArray->Reset();
178 while((candidate = static_cast<Candidate*>(fItTrackInputArray->Next())))
179 {
180 if (candidate->IsPU)
181 continue;
182 for (const auto &clusterID : clusterIDToIndex)
183 {
184 if (candidate->ClusterIndex != clusterID.first)
185 continue;
186 clusterIDToSumPT2.at (clusterID.first) += candidate->Momentum.Pt () * candidate->Momentum.Pt ();
187 }
188 }
189
190 for (const auto &clusterID : clusterIDToSumPT2)
191 sortedClusterIDs.push_back (make_pair (clusterID.first, clusterID.second));
192 sort (sortedClusterIDs.begin (), sortedClusterIDs.end (), secondDescending);
193 }
194 else
195 {
196 cout << "\"" << fMethod << "\" is not a valid sorting method!" << endl;
197 cout << "Valid methods are:" << endl;
198 cout << " BTV" << endl;
199 cout << " GenClosest" << endl;
200 cout << " GenBest" << endl;
201 throw 0;
202 }
203
204 for (const auto &clusterID : sortedClusterIDs)
205 {
206 Candidate *cluster = (Candidate *) fInputArray->At (clusterIDToIndex.at (clusterID.first));
207 if (fMethod == "BTV")
208 cluster->BTVSumPT2 = clusterID.second;
209 else if (fMethod == "GenClosest")
210 cluster->GenDeltaZ = clusterID.second;
211 else if (fMethod == "GenBest")
212 cluster->GenSumPT2 = clusterID.second;
213 fOutputArray->Add (cluster);
214 }
215}
216
217//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.