Fork me on GitHub

source: git/modules/VertexSorter.cc@ 5658083

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

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