Fork me on GitHub

source: git/modules/RunPUPPI.cc@ 4006893

Last change on this file since 4006893 was 341014c, checked in by Pavel Demin <pavel-demin@…>, 5 years ago

apply .clang-format to all .h, .cc and .cpp files
  • Property mode set to 100644
File size: 10.7 KB
Line 
1#include "modules/RunPUPPI.h"
2
3#include "PUPPI/AlgoObj.hh"
4#include "PUPPI/PuppiContainer.hh"
5#include "PUPPI/RecoObj2.hh"
6
7#include "fastjet/PseudoJet.hh"
8
9#include "classes/DelphesClasses.h"
10#include "classes/DelphesFactory.h"
11#include "classes/DelphesFormula.h"
12
13#include <algorithm>
14#include <iostream>
15#include <sstream>
16#include <stdexcept>
17#include <vector>
18
19using namespace std;
20using namespace fastjet;
21
22//------------------------------------------------------------------------------
23RunPUPPI::RunPUPPI() :
24 fItTrackInputArray(0),
25 fItNeutralInputArray(0)
26{
27}
28
29//------------------------------------------------------------------------------
30RunPUPPI::~RunPUPPI() {}
31
32//------------------------------------------------------------------------------
33
34void RunPUPPI::Init()
35{
36 // input collection
37 fTrackInputArray = ImportArray(GetString("TrackInputArray", "Calorimeter/towers"));
38 fItTrackInputArray = fTrackInputArray->MakeIterator();
39 fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "Calorimeter/towers"));
40 fItNeutralInputArray = fNeutralInputArray->MakeIterator();
41 fPVInputArray = ImportArray(GetString("PVInputArray", "PV"));
42 fPVItInputArray = fPVInputArray->MakeIterator();
43 // puppi parameters
44 fApplyNoLep = GetBool("UseNoLep", true);
45 fMinPuppiWeight = GetDouble("MinPuppiWeight", 0.01);
46 fUseExp = GetBool("UseExp", false);
47 // read eta min ranges
48 ExRootConfParam param = GetParam("EtaMinBin");
49 fEtaMinBin.clear();
50 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fEtaMinBin.push_back(param[iMap].GetDouble());
51 // read eta max ranges
52 param = GetParam("EtaMaxBin");
53 fEtaMaxBin.clear();
54 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fEtaMaxBin.push_back(param[iMap].GetDouble());
55 // read pt min value
56 param = GetParam("PtMinBin");
57 fPtMinBin.clear();
58 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fPtMinBin.push_back(param[iMap].GetDouble());
59 // read cone size
60 param = GetParam("ConeSizeBin");
61 fConeSizeBin.clear();
62 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fConeSizeBin.push_back(param[iMap].GetDouble());
63 // read RMS min pt
64 param = GetParam("RMSPtMinBin");
65 fRMSPtMinBin.clear();
66 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fRMSPtMinBin.push_back(param[iMap].GetDouble());
67 // read RMS scale factor
68 param = GetParam("RMSScaleFactorBin");
69 fRMSScaleFactorBin.clear();
70 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fRMSScaleFactorBin.push_back(param[iMap].GetDouble());
71 // read neutral pt min cut
72 param = GetParam("NeutralMinEBin");
73 fNeutralMinEBin.clear();
74 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralMinEBin.push_back(param[iMap].GetDouble());
75 // read neutral pt min slope
76 param = GetParam("NeutralPtSlope");
77 fNeutralPtSlope.clear();
78 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralPtSlope.push_back(param[iMap].GetDouble());
79 // read apply chs
80 //param = GetParam("ApplyCHS");
81 //fApplyCHS.clear();
82 //for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyCHS.push_back(param[iMap].GetBool());
83 // read use charged
84 param = GetParam("UseCharged");
85 fUseCharged.clear();
86 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fUseCharged.push_back(param[iMap].GetBool());
87 // read apply chs correction
88 param = GetParam("ApplyLowPUCorr");
89 fApplyLowPUCorr.clear();
90 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyLowPUCorr.push_back(param[iMap].GetBool());
91 // read metric id
92 param = GetParam("MetricId");
93 fMetricId.clear();
94 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fMetricId.push_back(param[iMap].GetInt());
95 // scheme for combining
96 param = GetParam("CombId");
97 fCombId.clear();
98 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fCombId.push_back(param[iMap].GetInt());
99 // create output array
100 fOutputArray = ExportArray(GetString("OutputArray", "puppiParticles"));
101 fOutputTrackArray = ExportArray(GetString("OutputArrayTracks", "puppiTracks"));
102 fOutputNeutralArray = ExportArray(GetString("OutputArrayNeutrals", "puppiNeutrals"));
103 // Create algorithm list for puppi
104 std::vector<AlgoObj> puppiAlgo;
105 if(puppiAlgo.empty())
106 {
107 if(!(fEtaMinBin.size() == fEtaMaxBin.size() and fEtaMinBin.size() == fPtMinBin.size() and fEtaMinBin.size() == fConeSizeBin.size() and fEtaMinBin.size() == fRMSPtMinBin.size()
108 and fEtaMinBin.size() == fRMSScaleFactorBin.size() and fEtaMinBin.size() == fNeutralMinEBin.size() and fEtaMinBin.size() == fNeutralPtSlope.size()
109 and fEtaMinBin.size() == fUseCharged.size()
110 and fEtaMinBin.size() == fApplyLowPUCorr.size() and fEtaMinBin.size() == fMetricId.size()))
111 {
112 std::cerr << " Error in PUPPI configuration, algo info should have the same size --> exit from the code" << std::endl;
113 std::exit(EXIT_FAILURE);
114 }
115 }
116 for(size_t iAlgo = 0; iAlgo < fEtaMinBin.size(); iAlgo++)
117 {
118 AlgoObj algoTmp;
119 algoTmp.etaMin = fEtaMinBin.at(iAlgo);
120 algoTmp.etaMax = fEtaMaxBin.at(iAlgo);
121 algoTmp.ptMin = fPtMinBin.at(iAlgo);
122 algoTmp.minNeutralPt = fNeutralMinEBin.at(iAlgo);
123 algoTmp.minNeutralPtSlope = fNeutralPtSlope.at(iAlgo);
124 //Eta Extrapolation stuff is missing
125 //Loop through file requiring algos for same bins to be adjacent
126 while(iAlgo < fEtaMinBin.size() and algoTmp.etaMin == fEtaMinBin.at(iAlgo) and algoTmp.etaMax == fEtaMaxBin.at(iAlgo))
127 {
128 AlgoSubObj algoSubTmp;
129 algoSubTmp.metricId = fMetricId.at(iAlgo);
130 algoSubTmp.useCharged = fUseCharged.at(iAlgo);
131 algoSubTmp.applyLowPUCorr = fApplyLowPUCorr.at(iAlgo);
132 algoSubTmp.combId = fCombId.at(iAlgo);
133 algoSubTmp.coneSize = fConeSizeBin.at(iAlgo);
134 algoSubTmp.rmsPtMin = fRMSPtMinBin.at(iAlgo);
135 algoSubTmp.rmsScaleFactor = fRMSScaleFactorBin.at(iAlgo);
136 algoTmp.subAlgos.push_back(algoSubTmp);
137 iAlgo++;
138 }
139 iAlgo--;
140 //if(std::find(puppiAlgo.begin(),puppiAlgo.end(),algoTmp) != puppiAlgo.end()) continue;
141 puppiAlgo.push_back(algoTmp);
142 }
143 fPuppi = new PuppiContainer(true, fUseExp, fMinPuppiWeight, puppiAlgo);
144}
145
146//------------------------------------------------------------------------------
147
148void RunPUPPI::Finish()
149{
150 if(fItTrackInputArray) delete fItTrackInputArray;
151 if(fItNeutralInputArray) delete fItNeutralInputArray;
152 if(fPuppi) delete fPuppi;
153}
154
155//------------------------------------------------------------------------------
156
157void RunPUPPI::Process()
158{
159
160 Candidate *candidate, *particle;
161 TLorentzVector momentum;
162
163 //DelphesFactory *factory = GetFactory();
164
165 // loop over input objects
166 fItTrackInputArray->Reset();
167 fItNeutralInputArray->Reset();
168 fPVItInputArray->Reset();
169
170 std::vector<Candidate *> InputParticles;
171 InputParticles.clear();
172
173 // take the leading vertex
174 float PVZ = 0.;
175 Candidate *pv = static_cast<Candidate *>(fPVItInputArray->Next());
176 if(pv) PVZ = pv->Position.Z();
177 // Fill input particles for puppi
178 std::vector<RecoObj> puppiInputVector;
179 puppiInputVector.clear();
180 int lNBad = 0;
181 // Loop on charge track candidate
182 while((candidate = static_cast<Candidate *>(fItTrackInputArray->Next())))
183 {
184 momentum = candidate->Momentum;
185 RecoObj curRecoObj;
186 curRecoObj.pt = momentum.Pt();
187 curRecoObj.eta = momentum.Eta();
188 curRecoObj.phi = momentum.Phi();
189 curRecoObj.m = momentum.M();
190 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0)); //if(fApplyNoLep && TMath::Abs(candidate->PID) == 11) continue; //Dumb cut to minimize the nolepton on electron
191 //if(fApplyNoLep && TMath::Abs(candidate->PID) == 13) continue;
192 if(candidate->IsRecoPU and candidate->Charge != 0)
193 { // if it comes fromPU vertexes after the resolution smearing and the dZ matching within resolution
194 lNBad++;
195 curRecoObj.id = 2;
196 curRecoObj.vtxId = 0.7 * (fPVInputArray->GetEntries()); //Hack apply reco vtx efficiency of 70% for calibration
197 if(TMath::Abs(candidate->PID) == 11)
198 curRecoObj.pfType = 2;
199 else if(TMath::Abs(candidate->PID) == 13)
200 curRecoObj.pfType = 3;
201 else if(TMath::Abs(candidate->PID) == 22)
202 curRecoObj.pfType = 4;
203 else
204 curRecoObj.pfType = 1;
205 curRecoObj.dZ = particle->Position.Z() - PVZ;
206 }
207 else if(!candidate->IsRecoPU && candidate->Charge != 0)
208 {
209 curRecoObj.id = 1; // charge from LV
210 curRecoObj.vtxId = 1; // from PV
211 if(TMath::Abs(candidate->PID) == 11)
212 curRecoObj.pfType = 2;
213 else if(TMath::Abs(candidate->PID) == 13)
214 curRecoObj.pfType = 3;
215 else if(TMath::Abs(candidate->PID) == 22)
216 curRecoObj.pfType = 4;
217 else
218 curRecoObj.pfType = 1;
219 curRecoObj.dZ = particle->Position.Z() - PVZ;
220 }
221 else
222 {
223 std::cerr << " RunPUPPI: problem with a charged track --> it has charge 0 " << std::endl;
224 continue;
225 }
226
227 puppiInputVector.push_back(curRecoObj);
228 InputParticles.push_back(candidate);
229 }
230
231 // Loop on neutral calo cells
232 while((candidate = static_cast<Candidate *>(fItNeutralInputArray->Next())))
233 {
234 momentum = candidate->Momentum;
235 RecoObj curRecoObj;
236 curRecoObj.pt = momentum.Pt();
237 curRecoObj.eta = momentum.Eta();
238 curRecoObj.phi = momentum.Phi();
239 curRecoObj.m = momentum.M();
240 curRecoObj.charge = 0;
241 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
242 if(candidate->Charge == 0)
243 {
244 curRecoObj.id = 0; // neutrals have id==0
245 curRecoObj.vtxId = 0; // neutrals have vtxId==0
246 if(TMath::Abs(candidate->PID) == 11)
247 curRecoObj.pfType = 2;
248 else if(TMath::Abs(candidate->PID) == 13)
249 curRecoObj.pfType = 3;
250 else if(TMath::Abs(candidate->PID) == 22)
251 curRecoObj.pfType = 4;
252 else
253 curRecoObj.pfType = 5;
254 curRecoObj.dZ = particle->Position.Z() - PVZ;
255 }
256 else
257 {
258 std::cerr << " RunPUPPI: problem with a neutrals cells --> it has charge !=0 " << std::endl;
259 continue;
260 }
261 puppiInputVector.push_back(curRecoObj);
262 InputParticles.push_back(candidate);
263 }
264 // Create PUPPI container
265 fPuppi->initialize(puppiInputVector);
266 fPuppi->puppiWeights();
267 std::vector<PseudoJet> puppiParticles = fPuppi->puppiParticles();
268
269 // Loop on final particles
270 for(std::vector<PseudoJet>::iterator it = puppiParticles.begin(); it != puppiParticles.end(); it++)
271 {
272 if(it->user_index() <= int(InputParticles.size()))
273 {
274 candidate = static_cast<Candidate *>(InputParticles.at(it->user_index())->Clone());
275 candidate->Momentum.SetPxPyPzE(it->px(), it->py(), it->pz(), it->e());
276 fOutputArray->Add(candidate);
277 if(puppiInputVector.at(it->user_index()).id == 1 or puppiInputVector.at(it->user_index()).id == 2)
278 fOutputTrackArray->Add(candidate);
279 else if(puppiInputVector.at(it->user_index()).id == 0)
280 fOutputNeutralArray->Add(candidate);
281 }
282 else
283 {
284 std::cerr << " particle not found in the input Array --> skip " << std::endl;
285 continue;
286 }
287 }
288}
Note: See TracBrowser for help on using the repository browser.