Fork me on GitHub

source: git/modules/RunPUPPI.cc@ 9da65a5

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

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