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 |
|
---|
19 | using namespace std;
|
---|
20 | using namespace fastjet;
|
---|
21 |
|
---|
22 | //------------------------------------------------------------------------------
|
---|
23 | RunPUPPI::RunPUPPI() :
|
---|
24 | fItTrackInputArray(0),
|
---|
25 | fItNeutralInputArray(0)
|
---|
26 | {}
|
---|
27 |
|
---|
28 | //------------------------------------------------------------------------------
|
---|
29 | RunPUPPI::~RunPUPPI(){}
|
---|
30 |
|
---|
31 | //------------------------------------------------------------------------------
|
---|
32 |
|
---|
33 | void 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 |
|
---|
142 | void RunPUPPI::Finish(){
|
---|
143 | if(fItTrackInputArray) delete fItTrackInputArray;
|
---|
144 | if(fItNeutralInputArray) delete fItNeutralInputArray;
|
---|
145 | }
|
---|
146 |
|
---|
147 | //------------------------------------------------------------------------------
|
---|
148 |
|
---|
149 | void 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 | }
|
---|