Fork me on GitHub

source: git/modules/TrackTimingPileUpSubtractor.cc@ b4a3c55

Timing
Last change on this file since b4a3c55 was 584e2e9, checked in by Kaan Yüksel Oyulmaz <kaanyukseloyulmaz@…>, 5 years ago

FCC-hh Vertex Card is updated with PileUpMerger

  • Property mode set to 100644
File size: 5.6 KB
Line 
1/*
2 * Delphes: a framework for fast simulation of a generic collider experiment
3 * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19/** \class TrackPileUpSubtractor
20 *
21 * Subtract pile-up contribution from tracks.
22 *
23 * \author P. Demin - UCL, Louvain-la-Neuve
24 *
25 */
26
27#include "modules/TrackTimingPileUpSubtractor.h"
28
29#include "classes/DelphesClasses.h"
30#include "classes/DelphesFactory.h"
31#include "classes/DelphesFormula.h"
32
33#include "ExRootAnalysis/ExRootClassifier.h"
34#include "ExRootAnalysis/ExRootFilter.h"
35#include "ExRootAnalysis/ExRootResult.h"
36
37#include "TDatabasePDG.h"
38#include "TFormula.h"
39#include "TLorentzVector.h"
40#include "TMath.h"
41#include "TObjArray.h"
42#include "TRandom3.h"
43#include "TString.h"
44
45#include <algorithm>
46#include <iostream>
47#include <sstream>
48#include <stdexcept>
49
50using namespace std;
51
52//------------------------------------------------------------------------------
53
54TrackTimingPileUpSubtractor::TrackTimingPileUpSubtractor() :
55 fFormula(0)
56{
57 fFormula = new DelphesFormula;
58}
59
60//------------------------------------------------------------------------------
61
62TrackTimingPileUpSubtractor::~TrackTimingPileUpSubtractor()
63{
64 if(fFormula) delete fFormula;
65}
66
67//------------------------------------------------------------------------------
68
69void TrackTimingPileUpSubtractor::Init()
70{
71 // import input array
72
73 fVertexInputArray = ImportArray(GetString("VertexInputArray", "PileUpMerger/vertices"));
74 fItVertexInputArray = fVertexInputArray->MakeIterator();
75
76 // read resolution formula in m
77 fFormula->Compile(GetString("ZVertexResolution", "0.001"));
78
79 fPTMin = GetDouble("PTMin", 0.);
80
81 // import arrays with output from other modules
82
83 ExRootConfParam param = GetParam("InputArray");
84 Long_t i, size;
85 const TObjArray *array;
86 TIterator *iterator;
87
88 size = param.GetSize();
89 for(i = 0; i < size / 2; ++i)
90 {
91 array = ImportArray(param[i * 2].GetString());
92 iterator = array->MakeIterator();
93
94 fInputMap[iterator] = ExportArray(param[i * 2 + 1].GetString());
95 }
96}
97
98//------------------------------------------------------------------------------
99
100void TrackTimingPileUpSubtractor::Finish()
101{
102 map<TIterator *, TObjArray *>::iterator itInputMap;
103 TIterator *iterator;
104
105 for(itInputMap = fInputMap.begin(); itInputMap != fInputMap.end(); ++itInputMap)
106 {
107 iterator = itInputMap->first;
108
109 if(iterator) delete iterator;
110 }
111
112 if(fItVertexInputArray) delete fItVertexInputArray;
113}
114
115//------------------------------------------------------------------------------
116
117void TrackTimingPileUpSubtractor::Process()
118{
119 Candidate *candidate, *particle;
120 map<TIterator *, TObjArray *>::iterator itInputMap;
121 TIterator *iterator;
122 TObjArray *array;
123 Double_t z, zvtx = 0;
124 Double_t z_err, zvtx_err = 0;
125 Double_t t, tvtx = 0;
126 Double_t t_err, tvtx_err = 0;
127 Double_t sumPTSquare = 0;
128 Double_t tempPTSquare = 0;
129 Double_t pt, eta, phi, e;
130 Double_t distance = 0;
131
132 // find z position of primary vertex
133
134 fItVertexInputArray->Reset();
135 while((candidate = static_cast<Candidate *>(fItVertexInputArray->Next())))
136 {
137 tempPTSquare = candidate->SumPT2;
138 if(tempPTSquare > sumPTSquare)
139 {
140 sumPTSquare = tempPTSquare;
141 zvtx = candidate->Position.Z();
142 zvtx_err = candidate->PositionError.Z();
143 tvtx = candidate->Position.T();
144 tvtx_err = candidate->PositionError.T();
145 cout << " initial : " << candidate->InitialPosition.T() << " final : " << candidate->Position.T() << endl;
146 }
147 }
148
149 // loop over all input arrays
150 for(itInputMap = fInputMap.begin(); itInputMap != fInputMap.end(); ++itInputMap)
151 {
152 iterator = itInputMap->first;
153 array = itInputMap->second;
154
155 // loop over all candidates
156 iterator->Reset();
157 while((candidate = static_cast<Candidate *>(iterator->Next())))
158 {
159 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
160 const TLorentzVector &candidateMomentum = particle->Momentum;
161
162 eta = candidateMomentum.Eta();
163 pt = candidateMomentum.Pt();
164 phi = candidateMomentum.Phi();
165 e = candidateMomentum.E();
166
167 z = particle->Position.Z();
168 z_err = particle->PositionError.Z();
169 t = particle->Position.T();
170 t_err = particle->PositionError.T();
171
172 // apply pile-up subtraction
173 distance = pow((zvtx - z),2)/pow((zvtx_err - z_err),2) + pow((tvtx - t),2)/pow((tvtx_err - t_err),2);
174 //cout << " t : " << tvtx << " t(particle)" << t << endl;
175 // here I calculated distance using Z and T of selected vertex (highest sum Pt square) and particles
176 // however z_err of vertices is gives 0 because of using CMS trackResolutionCMS.tcl (in that formula, there is limitation on |eta| < 2.5)
177 // thats why I used TMath::Abs(z - zvtx) < 0.005 && TMath::Abs(t - tvtx) < 5.0
178
179 if(candidate->Charge != 0 && TMath::Abs(z - zvtx) < 0.005 && TMath::Abs(t - tvtx) < 5.0)
180 {
181 candidate->IsRecoPU = 1;
182 }
183 else
184 {
185 candidate->IsRecoPU = 0;
186 if(candidate->Momentum.Pt() > fPTMin) array->Add(candidate);
187 }
188 }
189 }
190}
Note: See TracBrowser for help on using the repository browser.