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