- Timestamp:
- Jun 6, 2016, 12:16:53 PM (8 years ago)
- Branches:
- ImprovedOutputFile, Timing, dual_readout, llp, master
- Children:
- 94c5375
- Parents:
- 989ee8f
- Location:
- modules
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/RunPUPPI.cc
r989ee8f re921a28 1 1 #include "modules/RunPUPPI.h" 2 2 3 #include "PUPPI/ puppiCleanContainer.hh"4 #include "PUPPI/ RecoObj.hh"5 #include "PUPPI/puppiParticle.hh"6 #include "PUPPI/puppiAlgoBin.hh"3 #include "PUPPI/RecoObj2.hh" 4 #include "PUPPI/AlgoObj.hh" 5 //#include "PUPPI/puppiParticle.hh" 6 //#include "PUPPI/puppiAlgoBin.hh" 7 7 8 8 #include "classes/DelphesClasses.h" … … 30 30 31 31 void RunPUPPI::Init(){ 32 33 32 // input collection 34 33 fTrackInputArray = ImportArray(GetString("TrackInputArray", "Calorimeter/towers")); … … 38 37 fPVInputArray = ImportArray(GetString("PVInputArray", "PV")); 39 38 fPVItInputArray = fPVInputArray->MakeIterator(); 40 41 42 // puppi parameters 39 // puppi parameters 40 fApplyNoLep = GetBool("UseNoLep", true); 43 41 fMinPuppiWeight = GetDouble("MinPuppiWeight", 0.01); 44 42 fUseExp = GetBool("UseExp", false); 45 46 // read eta min ranges 43 // read eta min ranges 47 44 ExRootConfParam param = GetParam("EtaMinBin"); 48 45 fEtaMinBin.clear(); 49 46 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fEtaMinBin.push_back(param[iMap].GetDouble()); 50 51 // read eta max ranges 47 // read eta max ranges 52 48 param = GetParam("EtaMaxBin"); 53 49 fEtaMaxBin.clear(); 54 50 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fEtaMaxBin.push_back(param[iMap].GetDouble()); 55 56 // read pt min value 51 // read pt min value 57 52 param = GetParam("PtMinBin"); 58 53 fPtMinBin.clear(); 59 54 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fPtMinBin.push_back(param[iMap].GetDouble()); 60 61 55 // read cone size 62 56 param = GetParam("ConeSizeBin"); 63 57 fConeSizeBin.clear(); 64 58 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fConeSizeBin.push_back(param[iMap].GetDouble()); 65 66 59 // read RMS min pt 67 60 param = GetParam("RMSPtMinBin"); 68 61 fRMSPtMinBin.clear(); 69 62 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fRMSPtMinBin.push_back(param[iMap].GetDouble()); 70 71 // read RMS scale factor 63 // read RMS scale factor 72 64 param = GetParam("RMSScaleFactorBin"); 73 65 fRMSScaleFactorBin.clear(); 74 66 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fRMSScaleFactorBin.push_back(param[iMap].GetDouble()); 75 76 // read neutral pt min cut 67 // read neutral pt min cut 77 68 param = GetParam("NeutralMinEBin"); 78 69 fNeutralMinEBin.clear(); 79 70 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralMinEBin.push_back(param[iMap].GetDouble()); 80 81 // read neutral pt min slope 71 // read neutral pt min slope 82 72 param = GetParam("NeutralPtSlope"); 83 73 fNeutralPtSlope.clear(); 84 74 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralPtSlope.push_back(param[iMap].GetDouble()); 85 86 75 // 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 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 92 80 param = GetParam("UseCharged"); 93 81 fUseCharged.clear(); 94 82 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fUseCharged.push_back(param[iMap].GetBool()); 95 96 // read apply chs correction 83 // read apply chs correction 97 84 param = GetParam("ApplyLowPUCorr"); 98 85 fApplyLowPUCorr.clear(); 99 86 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyLowPUCorr.push_back(param[iMap].GetBool()); 100 101 87 // read metric id 102 88 param = GetParam("MetricId"); 103 89 fMetricId.clear(); 104 90 for(int iMap = 0; iMap < param.GetSize(); ++iMap) fMetricId.push_back(param[iMap].GetInt()); 105 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()); 106 95 // create output array 107 96 fOutputArray = ExportArray(GetString("OutputArray", "puppiParticles")); 108 97 fOutputTrackArray = ExportArray(GetString("OutputArrayTracks", "puppiTracks")); 109 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); 110 136 } 111 137 … … 124 150 TLorentzVector momentum; 125 151 126 DelphesFactory *factory = GetFactory();152 //DelphesFactory *factory = GetFactory(); 127 153 128 154 // loop over input objects 129 fItTrackInputArray ->Reset();130 fItNeutralInputArray ->Reset();131 fPVItInputArray ->Reset();155 fItTrackInputArray ->Reset(); 156 fItNeutralInputArray ->Reset(); 157 fPVItInputArray ->Reset(); 132 158 133 159 std::vector<Candidate *> InputParticles; … … 142 168 std::vector<RecoObj> puppiInputVector; 143 169 puppiInputVector.clear(); 144 170 int lNBad = 0; 145 171 // Loop on charge track candidate 146 172 while((candidate = static_cast<Candidate*>(fItTrackInputArray->Next()))){ 147 148 173 momentum = candidate->Momentum; 149 150 174 RecoObj curRecoObj; 151 175 curRecoObj.pt = momentum.Pt(); … … 154 178 curRecoObj.m = momentum.M(); 155 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; 156 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++; 157 184 curRecoObj.id = 2; 158 curRecoObj.vtxId = candidate->IsPU;185 curRecoObj.vtxId = 0.7*(fPVInputArray->GetEntries()); //Hack apply reco vtx efficiency of 70% for calibration 159 186 if(TMath::Abs(candidate->PID) == 11) curRecoObj.pfType = 2; 160 187 else if(TMath::Abs(candidate->PID) == 13) curRecoObj.pfType = 3; … … 183 210 // Loop on neutral calo cells 184 211 while((candidate = static_cast<Candidate*>(fItNeutralInputArray->Next()))){ 185 186 212 momentum = candidate->Momentum; 187 188 213 RecoObj curRecoObj; 189 214 curRecoObj.pt = momentum.Pt(); … … 191 216 curRecoObj.phi = momentum.Phi(); 192 217 curRecoObj.m = momentum.M(); 218 curRecoObj.charge = 0; 193 219 particle = static_cast<Candidate*>(candidate->GetCandidates()->Last()); 194 195 220 196 221 if(candidate->Charge == 0){ … … 210 235 InputParticles.push_back(candidate); 211 236 } 212 213 // Create algorithm list for puppi214 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 237 // Create PUPPI container 244 puppiCleanContainer curEvent(puppiInputVector,puppiAlgo,fMinPuppiWeight,fUseExp); 245 std::vector<fastjet::PseudoJet> puppiParticles = curEvent.puppiEvent(); 238 fPuppi->initialize(puppiInputVector); 239 fPuppi->puppiWeights(); 240 std::vector<fastjet::PseudoJet> puppiParticles = fPuppi->puppiParticles(); 246 241 247 242 // Loop on final particles -
modules/RunPUPPI.h
r989ee8f re921a28 3 3 4 4 #include "classes/DelphesModule.h" 5 #include "PUPPI/PuppiContainer.hh" 5 6 #include <vector> 6 7 … … 29 30 const TObjArray *fNeutralInputArray; //! 30 31 const TObjArray *fPVInputArray; //! 31 32 PuppiContainer* fPuppi; 32 33 // puppi parameters 33 float fMinPuppiWeight; 34 bool fApplyNoLep; 35 double fMinPuppiWeight; 34 36 bool fUseExp; 35 37 … … 46 48 std::vector<bool> fApplyLowPUCorr; 47 49 std::vector<int> fMetricId; 50 std::vector<int> fCombId; 48 51 49 52 TObjArray *fOutputArray;
Note:
See TracChangeset
for help on using the changeset viewer.