Fork me on GitHub

source: git/modules/VertexSorter.cc @ 77e9ae1

ImprovedOutputFileTimingllp
Last change on this file since 77e9ae1 was 77e9ae1, checked in by Pavel Demin <pavel-demin@…>, 18 months 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.