Fork me on GitHub

source: git/modules/RunPUPPI.cc@ a393b57

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

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