Fork me on GitHub

source: git/modules/RunPUPPI.cc@ ea0f50d

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

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