Changeset 341014c in git for modules/PileUpJetID.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/PileUpJetID.cc
r45e58be r341014c 14 14 #include "classes/DelphesFormula.h" 15 15 16 #include "ExRootAnalysis/ExRootClassifier.h" 17 #include "ExRootAnalysis/ExRootFilter.h" 16 18 #include "ExRootAnalysis/ExRootResult.h" 17 #include "ExRootAnalysis/ExRootFilter.h" 18 #include "ExRootAnalysis/ExRootClassifier.h" 19 19 20 #include "TFormula.h" 20 21 #include "TMath.h" 22 #include "TObjArray.h" 23 #include "TRandom3.h" 21 24 #include "TString.h" 22 #include "TFormula.h"23 #include "TRandom3.h"24 #include "TObjArray.h"25 25 //#include "TDatabasePDG.h" 26 26 #include "TLorentzVector.h" 27 27 28 28 #include <algorithm> 29 #include <stdexcept>30 29 #include <iostream> 31 30 #include <sstream> 31 #include <stdexcept> 32 32 33 33 using namespace std; … … 36 36 37 37 PileUpJetID::PileUpJetID() : 38 fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0) 39 { 40 38 fItJetInputArray(0), fTrackInputArray(0), fNeutralInputArray(0) 39 { 41 40 } 42 41 … … 45 44 PileUpJetID::~PileUpJetID() 46 45 { 47 48 46 } 49 47 … … 57 55 fUseConstituents = GetInt("UseConstituents", 0); 58 56 59 fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel", 0.1);60 fBetaMinBarrel = GetDouble("BetaMinBarrel", 0.1);61 fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap", 0.1);62 fBetaMinEndcap = GetDouble("BetaMinEndcap", 0.1);63 fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward", 0.1);57 fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel", 0.1); 58 fBetaMinBarrel = GetDouble("BetaMinBarrel", 0.1); 59 fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap", 0.1); 60 fBetaMinEndcap = GetDouble("BetaMinEndcap", 0.1); 61 fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward", 0.1); 64 62 fJetPTMinForNeutrals = GetDouble("JetPTMinForNeutrals", 20.0); 65 63 fNeutralPTMin = GetDouble("NeutralPTMin", 2.0); … … 82 80 fOutputArray = ExportArray(GetString("OutputArray", "jets")); 83 81 84 fNeutralsInPassingJets = ExportArray(GetString("NeutralsInPassingJets","eflowtowers")); 85 82 fNeutralsInPassingJets = ExportArray(GetString("NeutralsInPassingJets", "eflowtowers")); 86 83 } 87 84 … … 95 92 if(fItTrackInputArray) delete fItTrackInputArray; 96 93 if(fItNeutralInputArray) delete fItNeutralInputArray; 97 98 94 } 99 95 … … 109 105 // loop over all input candidates 110 106 fItJetInputArray->Reset(); 111 while((candidate = static_cast<Candidate *>(fItJetInputArray->Next())))107 while((candidate = static_cast<Candidate *>(fItJetInputArray->Next()))) 112 108 { 113 109 momentum = candidate->Momentum; … … 133 129 float pt_ann[5]; 134 130 135 for (int i = 0 ; i < 5 ; i++) { 131 for(int i = 0; i < 5; i++) 132 { 136 133 pt_ann[i] = 0.; 137 134 } 138 135 139 if (fUseConstituents) { 136 if(fUseConstituents) 137 { 140 138 TIter itConstituents(candidate->GetCandidates()); 141 while((constituent = static_cast<Candidate*>(itConstituents.Next()))) { 139 while((constituent = static_cast<Candidate *>(itConstituents.Next()))) 140 { 142 141 float pt = constituent->Momentum.Pt(); 143 142 float dr = candidate->Momentum.DeltaR(constituent->Momentum); 144 // cout << " There exists a constituent with dr=" << dr << endl; 145 sumpt += pt; 146 sumdrsqptsq += dr*dr*pt*pt; 147 sumptsq += pt*pt; 148 if (constituent->Charge == 0) { 149 nn++; 150 } else { 151 if (constituent->IsRecoPU) { 152 sumptchpu += pt; 153 } else { 154 sumptchpv += pt; 155 } 156 sumptch += pt; 157 nc++; 158 } 159 for (int i = 0 ; i < 5 ; i++) { 160 if (dr > 0.1*i && dr < 0.1*(i+1)) { 161 pt_ann[i] += pt; 162 } 163 } 164 float tow_sumT = 0; 165 float tow_sumW = 0; 166 for (int i = 0 ; i < constituent->ECalEnergyTimePairs.size() ; i++) { 167 float w = TMath::Sqrt(constituent->ECalEnergyTimePairs[i].first); 168 if (fAverageEachTower) { 169 tow_sumT += w*constituent->ECalEnergyTimePairs[i].second; 143 // cout << " There exists a constituent with dr=" << dr << endl; 144 sumpt += pt; 145 sumdrsqptsq += dr * dr * pt * pt; 146 sumptsq += pt * pt; 147 if(constituent->Charge == 0) 148 { 149 nn++; 150 } 151 else 152 { 153 if(constituent->IsRecoPU) 154 { 155 sumptchpu += pt; 156 } 157 else 158 { 159 sumptchpv += pt; 160 } 161 sumptch += pt; 162 nc++; 163 } 164 for(int i = 0; i < 5; i++) 165 { 166 if(dr > 0.1 * i && dr < 0.1 * (i + 1)) 167 { 168 pt_ann[i] += pt; 169 } 170 } 171 float tow_sumT = 0; 172 float tow_sumW = 0; 173 for(int i = 0; i < constituent->ECalEnergyTimePairs.size(); i++) 174 { 175 float w = TMath::Sqrt(constituent->ECalEnergyTimePairs[i].first); 176 if(fAverageEachTower) 177 { 178 tow_sumT += w * constituent->ECalEnergyTimePairs[i].second; 170 179 tow_sumW += w; 171 } else { 172 sumT0 += w*constituent->ECalEnergyTimePairs[i].second; 173 sumT1 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.001); 174 sumT10 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.010); 175 sumT20 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.020); 176 sumT30 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.030); 177 sumT40 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.040); 178 sumWeightsForT += w; 179 candidate->NTimeHits++; 180 } 181 } 182 if (fAverageEachTower && tow_sumW > 0.) { 183 sumT0 += tow_sumT; 184 sumT1 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.001); 185 sumT10 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0010); 186 sumT20 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0020); 187 sumT30 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0030); 188 sumT40 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0040); 189 sumWeightsForT += tow_sumW; 190 candidate->NTimeHits++; 191 } 192 } 193 } else { 180 } 181 else 182 { 183 sumT0 += w * constituent->ECalEnergyTimePairs[i].second; 184 sumT1 += w * gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second, 0.001); 185 sumT10 += w * gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second, 0.010); 186 sumT20 += w * gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second, 0.020); 187 sumT30 += w * gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second, 0.030); 188 sumT40 += w * gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second, 0.040); 189 sumWeightsForT += w; 190 candidate->NTimeHits++; 191 } 192 } 193 if(fAverageEachTower && tow_sumW > 0.) 194 { 195 sumT0 += tow_sumT; 196 sumT1 += tow_sumW * gRandom->Gaus(tow_sumT / tow_sumW, 0.001); 197 sumT10 += tow_sumW * gRandom->Gaus(tow_sumT / tow_sumW, 0.0010); 198 sumT20 += tow_sumW * gRandom->Gaus(tow_sumT / tow_sumW, 0.0020); 199 sumT30 += tow_sumW * gRandom->Gaus(tow_sumT / tow_sumW, 0.0030); 200 sumT40 += tow_sumW * gRandom->Gaus(tow_sumT / tow_sumW, 0.0040); 201 sumWeightsForT += tow_sumW; 202 candidate->NTimeHits++; 203 } 204 } 205 } 206 else 207 { 194 208 // Not using constituents, using dr 195 209 fItTrackInputArray->Reset(); 196 while ((trk = static_cast<Candidate*>(fItTrackInputArray->Next()))) { 197 if (trk->Momentum.DeltaR(candidate->Momentum) < fParameterR) { 198 float pt = trk->Momentum.Pt(); 199 sumpt += pt; 200 sumptch += pt; 201 if (trk->IsRecoPU) { 202 sumptchpu += pt; 203 } else { 204 sumptchpv += pt; 205 } 206 float dr = candidate->Momentum.DeltaR(trk->Momentum); 207 sumdrsqptsq += dr*dr*pt*pt; 208 sumptsq += pt*pt; 209 nc++; 210 for (int i = 0 ; i < 5 ; i++) { 211 if (dr > 0.1*i && dr < 0.1*(i+1)) { 212 pt_ann[i] += pt; 213 } 214 } 215 } 210 while((trk = static_cast<Candidate *>(fItTrackInputArray->Next()))) 211 { 212 if(trk->Momentum.DeltaR(candidate->Momentum) < fParameterR) 213 { 214 float pt = trk->Momentum.Pt(); 215 sumpt += pt; 216 sumptch += pt; 217 if(trk->IsRecoPU) 218 { 219 sumptchpu += pt; 220 } 221 else 222 { 223 sumptchpv += pt; 224 } 225 float dr = candidate->Momentum.DeltaR(trk->Momentum); 226 sumdrsqptsq += dr * dr * pt * pt; 227 sumptsq += pt * pt; 228 nc++; 229 for(int i = 0; i < 5; i++) 230 { 231 if(dr > 0.1 * i && dr < 0.1 * (i + 1)) 232 { 233 pt_ann[i] += pt; 234 } 235 } 236 } 216 237 } 217 238 fItNeutralInputArray->Reset(); 218 while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) { 219 if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR) { 220 float pt = constituent->Momentum.Pt(); 221 sumpt += pt; 222 float dr = candidate->Momentum.DeltaR(constituent->Momentum); 223 sumdrsqptsq += dr*dr*pt*pt; 224 sumptsq += pt*pt; 225 nn++; 226 for (int i = 0 ; i < 5 ; i++) { 227 if (dr > 0.1*i && dr < 0.1*(i+1)) { 228 pt_ann[i] += pt; 229 } 230 } 231 } 232 } 233 } 234 235 if (sumptch > 0.) { 236 candidate->Beta = sumptchpv/sumptch; 237 candidate->BetaStar = sumptchpu/sumptch; 238 } else { 239 while((constituent = static_cast<Candidate *>(fItNeutralInputArray->Next()))) 240 { 241 if(constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR) 242 { 243 float pt = constituent->Momentum.Pt(); 244 sumpt += pt; 245 float dr = candidate->Momentum.DeltaR(constituent->Momentum); 246 sumdrsqptsq += dr * dr * pt * pt; 247 sumptsq += pt * pt; 248 nn++; 249 for(int i = 0; i < 5; i++) 250 { 251 if(dr > 0.1 * i && dr < 0.1 * (i + 1)) 252 { 253 pt_ann[i] += pt; 254 } 255 } 256 } 257 } 258 } 259 260 if(sumptch > 0.) 261 { 262 candidate->Beta = sumptchpv / sumptch; 263 candidate->BetaStar = sumptchpu / sumptch; 264 } 265 else 266 { 239 267 candidate->Beta = -999.; 240 268 candidate->BetaStar = -999.; 241 269 } 242 if (sumptsq > 0.) { 243 candidate->MeanSqDeltaR = sumdrsqptsq/sumptsq; 244 } else { 270 if(sumptsq > 0.) 271 { 272 candidate->MeanSqDeltaR = sumdrsqptsq / sumptsq; 273 } 274 else 275 { 245 276 candidate->MeanSqDeltaR = -999.; 246 277 } 247 278 candidate->NCharged = nc; 248 279 candidate->NNeutrals = nn; 249 if (sumpt > 0.) { 280 if(sumpt > 0.) 281 { 250 282 candidate->PTD = TMath::Sqrt(sumptsq) / sumpt; 251 for (int i = 0 ; i < 5 ; i++) { 252 candidate->FracPt[i] = pt_ann[i]/sumpt; 253 } 254 } else { 283 for(int i = 0; i < 5; i++) 284 { 285 candidate->FracPt[i] = pt_ann[i] / sumpt; 286 } 287 } 288 else 289 { 255 290 candidate->PTD = -999.; 256 for (int i = 0 ; i < 5 ; i++) { 291 for(int i = 0; i < 5; i++) 292 { 257 293 candidate->FracPt[i] = -999.; 258 294 } … … 271 307 272 308 bool passId = false; 273 if (candidate->Momentum.Pt() > fJetPTMinForNeutrals && candidate->MeanSqDeltaR > -0.1) { 274 if (fabs(candidate->Momentum.Eta())<1.5) { 275 passId = ((candidate->Beta > fBetaMinBarrel) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxBarrel)); 276 } else if (fabs(candidate->Momentum.Eta())<4.0) { 277 passId = ((candidate->Beta > fBetaMinEndcap) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxEndcap)); 278 } else { 279 passId = (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxForward); 280 } 281 } 282 283 // cout << " Pt Eta MeanSqDeltaR Beta PassId " << candidate->Momentum.Pt() 309 if(candidate->Momentum.Pt() > fJetPTMinForNeutrals && candidate->MeanSqDeltaR > -0.1) 310 { 311 if(fabs(candidate->Momentum.Eta()) < 1.5) 312 { 313 passId = ((candidate->Beta > fBetaMinBarrel) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxBarrel)); 314 } 315 else if(fabs(candidate->Momentum.Eta()) < 4.0) 316 { 317 passId = ((candidate->Beta > fBetaMinEndcap) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxEndcap)); 318 } 319 else 320 { 321 passId = (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxForward); 322 } 323 } 324 325 // cout << " Pt Eta MeanSqDeltaR Beta PassId " << candidate->Momentum.Pt() 284 326 // << " " << candidate->Momentum.Eta() << " " << candidate->MeanSqDeltaR << " " << candidate->Beta << " " << passId << endl; 285 327 286 if (passId) { 287 if (fUseConstituents) { 288 TIter itConstituents(candidate->GetCandidates()); 289 while((constituent = static_cast<Candidate*>(itConstituents.Next()))) { 290 if (constituent->Charge == 0 && constituent->Momentum.Pt() > fNeutralPTMin) { 291 fNeutralsInPassingJets->Add(constituent); 292 // cout << " Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl; 293 } 294 } 295 } else { // use DeltaR 296 fItNeutralInputArray->Reset(); 297 while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) { 298 if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR && constituent->Momentum.Pt() > fNeutralPTMin) { 299 fNeutralsInPassingJets->Add(constituent); 300 // cout << " Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl; 301 } 302 } 303 } 304 } 305 306 328 if(passId) 329 { 330 if(fUseConstituents) 331 { 332 TIter itConstituents(candidate->GetCandidates()); 333 while((constituent = static_cast<Candidate *>(itConstituents.Next()))) 334 { 335 if(constituent->Charge == 0 && constituent->Momentum.Pt() > fNeutralPTMin) 336 { 337 fNeutralsInPassingJets->Add(constituent); 338 // cout << " Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl; 339 } 340 } 341 } 342 else 343 { // use DeltaR 344 fItNeutralInputArray->Reset(); 345 while((constituent = static_cast<Candidate *>(fItNeutralInputArray->Next()))) 346 { 347 if(constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR && constituent->Momentum.Pt() > fNeutralPTMin) 348 { 349 fNeutralsInPassingJets->Add(constituent); 350 // cout << " Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl; 351 } 352 } 353 } 354 } 307 355 } 308 356 }
Note:
See TracChangeset
for help on using the changeset viewer.