Changeset 341014c in git for modules/RunPUPPI.cc
- Timestamp:
- Feb 12, 2019, 9:29:17 PM (6 years ago)
- Branches:
- ImprovedOutputFile, Timing, llp, master
- Children:
- 6455202
- Parents:
- 45e58be
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/RunPUPPI.cc
r45e58be r341014c 1 1 #include "modules/RunPUPPI.h" 2 2 3 #include "PUPPI/RecoObj2.hh"4 3 #include "PUPPI/AlgoObj.hh" 5 4 #include "PUPPI/PuppiContainer.hh" 5 #include "PUPPI/RecoObj2.hh" 6 6 7 7 #include "fastjet/PseudoJet.hh" … … 11 11 #include "classes/DelphesFormula.h" 12 12 13 #include <algorithm> 14 #include <stdexcept> 13 #include <algorithm> 15 14 #include <iostream> 16 15 #include <sstream> 16 #include <stdexcept> 17 17 #include <vector> 18 18 … … 22 22 //------------------------------------------------------------------------------ 23 23 RunPUPPI::RunPUPPI() : 24 fItTrackInputArray(0), 24 fItTrackInputArray(0), 25 25 fItNeutralInputArray(0) 26 {} 27 28 //------------------------------------------------------------------------------ 29 RunPUPPI::~RunPUPPI(){} 30 31 //------------------------------------------------------------------------------ 32 33 void RunPUPPI::Init(){ 26 { 27 } 28 29 //------------------------------------------------------------------------------ 30 RunPUPPI::~RunPUPPI() {} 31 32 //------------------------------------------------------------------------------ 33 34 void RunPUPPI::Init() 35 { 34 36 // input collection 35 fTrackInputArray 36 fItTrackInputArray 37 fNeutralInputArray 37 fTrackInputArray = ImportArray(GetString("TrackInputArray", "Calorimeter/towers")); 38 fItTrackInputArray = fTrackInputArray->MakeIterator(); 39 fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "Calorimeter/towers")); 38 40 fItNeutralInputArray = fNeutralInputArray->MakeIterator(); 39 fPVInputArray 40 fPVItInputArray 41 // puppi parameters 42 fApplyNoLep = GetBool("UseNoLep", true);41 fPVInputArray = ImportArray(GetString("PVInputArray", "PV")); 42 fPVItInputArray = fPVInputArray->MakeIterator(); 43 // puppi parameters 44 fApplyNoLep = GetBool("UseNoLep", true); 43 45 fMinPuppiWeight = GetDouble("MinPuppiWeight", 0.01); 44 fUseExp 46 fUseExp = GetBool("UseExp", false); 45 47 // read eta min ranges 46 48 ExRootConfParam param = GetParam("EtaMinBin"); … … 55 57 fPtMinBin.clear(); 56 58 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fPtMinBin.push_back(param[iMap].GetDouble()); 57 // read cone size 59 // read cone size 58 60 param = GetParam("ConeSizeBin"); 59 61 fConeSizeBin.clear(); 60 62 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fConeSizeBin.push_back(param[iMap].GetDouble()); 61 // read RMS min pt 63 // read RMS min pt 62 64 param = GetParam("RMSPtMinBin"); 63 65 fRMSPtMinBin.clear(); … … 75 77 fNeutralPtSlope.clear(); 76 78 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralPtSlope.push_back(param[iMap].GetDouble()); 77 // read apply chs 79 // read apply chs 78 80 //param = GetParam("ApplyCHS"); 79 81 //fApplyCHS.clear(); … … 87 89 fApplyLowPUCorr.clear(); 88 90 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyLowPUCorr.push_back(param[iMap].GetBool()); 89 // read metric id 91 // read metric id 90 92 param = GetParam("MetricId"); 91 93 fMetricId.clear(); … … 96 98 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fCombId.push_back(param[iMap].GetInt()); 97 99 // create output array 98 fOutputArray 99 fOutputTrackArray 100 fOutputArray = ExportArray(GetString("OutputArray", "puppiParticles")); 101 fOutputTrackArray = ExportArray(GetString("OutputArrayTracks", "puppiTracks")); 100 102 fOutputNeutralArray = ExportArray(GetString("OutputArrayNeutrals", "puppiNeutrals")); 101 103 // Create algorithm list for puppi 102 104 std::vector<AlgoObj> puppiAlgo; 103 if(puppiAlgo.empty()){ 105 if(puppiAlgo.empty()) 106 { 104 107 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; 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; 109 113 std::exit(EXIT_FAILURE); 110 114 } 111 115 } 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); 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); 118 123 algoTmp.minNeutralPtSlope = fNeutralPtSlope.at(iAlgo); 119 124 //Eta Extrapolation stuff is missing 120 125 //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)) { 126 while(iAlgo < fEtaMinBin.size() and algoTmp.etaMin == fEtaMinBin.at(iAlgo) and algoTmp.etaMax == fEtaMaxBin.at(iAlgo)) 127 { 122 128 AlgoSubObj algoSubTmp; 123 algoSubTmp.metricId 124 algoSubTmp.useCharged 125 algoSubTmp.applyLowPUCorr 126 algoSubTmp.combId 127 algoSubTmp.coneSize 128 algoSubTmp.rmsPtMin 129 algoSubTmp.rmsScaleFactor 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); 130 136 algoTmp.subAlgos.push_back(algoSubTmp); 131 137 iAlgo++; 132 138 } 133 139 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; 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 148 void RunPUPPI::Finish() 149 { 150 if(fItTrackInputArray) delete fItTrackInputArray; 144 151 if(fItNeutralInputArray) delete fItNeutralInputArray; 145 152 if(fPuppi) delete fPuppi; … … 148 155 //------------------------------------------------------------------------------ 149 156 150 void RunPUPPI::Process(){ 157 void RunPUPPI::Process() 158 { 151 159 152 160 Candidate *candidate, *particle; … … 156 164 157 165 // loop over input objects 158 fItTrackInputArray 159 fItNeutralInputArray 160 fPVItInputArray 166 fItTrackInputArray->Reset(); 167 fItNeutralInputArray->Reset(); 168 fPVItInputArray->Reset(); 161 169 162 170 std::vector<Candidate *> InputParticles; 163 171 InputParticles.clear(); 164 172 165 // take the leading vertex 173 // take the leading vertex 166 174 float PVZ = 0.; 167 Candidate *pv = static_cast<Candidate *>(fPVItInputArray->Next());168 if 175 Candidate *pv = static_cast<Candidate *>(fPVItInputArray->Next()); 176 if(pv) PVZ = pv->Position.Z(); 169 177 // Fill input particles for puppi 170 178 std::vector<RecoObj> puppiInputVector; 171 179 puppiInputVector.clear(); 172 int lNBad = 0;180 int lNBad = 0; 173 181 // 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()->At(0));//if(fApplyNoLep && TMath::Abs(candidate->PID) == 11) continue; //Dumb cut to minimize the nolepton on electron 182 //if(fApplyNoLep && TMath::Abs(candidate->PID) == 13) continue; 183 if (candidate->IsRecoPU and candidate->Charge !=0) { // if it comes fromPU vertexes after the resolution smearing and the dZ matching within resolution 184 lNBad++; 185 curRecoObj.id = 2; 186 curRecoObj.vtxId = 0.7*(fPVInputArray->GetEntries()); //Hack apply reco vtx efficiency of 70% for calibration 187 if(TMath::Abs(candidate->PID) == 11) curRecoObj.pfType = 2; 188 else if(TMath::Abs(candidate->PID) == 13) curRecoObj.pfType = 3; 189 else if(TMath::Abs(candidate->PID) == 22) curRecoObj.pfType = 4; 190 else curRecoObj.pfType = 1; 191 curRecoObj.dZ = particle->Position.Z()-PVZ; 192 } 193 else if(!candidate->IsRecoPU && candidate->Charge !=0) { 194 curRecoObj.id = 1; // charge from LV 195 curRecoObj.vtxId = 1; // from PV 196 if(TMath::Abs(candidate->PID) == 11) curRecoObj.pfType = 2; 197 else if(TMath::Abs(candidate->PID) == 13) curRecoObj.pfType = 3; 198 else if(TMath::Abs(candidate->PID) == 22) curRecoObj.pfType = 4; 199 else curRecoObj.pfType = 1; 200 curRecoObj.dZ = particle->Position.Z()-PVZ; 201 } 202 else { 203 std::cerr<<" RunPUPPI: problem with a charged track --> it has charge 0 "<<std::endl; 204 continue; 205 } 206 207 puppiInputVector.push_back(curRecoObj); 208 InputParticles.push_back(candidate); 209 } 210 211 // Loop on neutral calo cells 212 while((candidate = static_cast<Candidate*>(fItNeutralInputArray->Next()))){ 213 momentum = candidate->Momentum; 214 RecoObj curRecoObj; 215 curRecoObj.pt = momentum.Pt(); 216 curRecoObj.eta = momentum.Eta(); 217 curRecoObj.phi = momentum.Phi(); 218 curRecoObj.m = momentum.M(); 219 curRecoObj.charge = 0; 220 particle = static_cast<Candidate*>(candidate->GetCandidates()->At(0)); 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); 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); 236 263 } 237 264 // Create PUPPI container … … 241 268 242 269 // Loop on final particles 243 for (std::vector<PseudoJet>::iterator it = puppiParticles.begin() ; it != puppiParticles.end() ; it++) { 244 if(it->user_index() <= int(InputParticles.size())){ 270 for(std::vector<PseudoJet>::iterator it = puppiParticles.begin(); it != puppiParticles.end(); it++) 271 { 272 if(it->user_index() <= int(InputParticles.size())) 273 { 245 274 candidate = static_cast<Candidate *>(InputParticles.at(it->user_index())->Clone()); 246 candidate->Momentum.SetPxPyPzE(it->px(), it->py(),it->pz(),it->e());275 candidate->Momentum.SetPxPyPzE(it->px(), it->py(), it->pz(), it->e()); 247 276 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; 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; 253 285 continue; 254 } 255 } 256 257 } 286 } 287 } 288 }
Note:
See TracChangeset
for help on using the changeset viewer.