File indexing completed on 2023-10-25 10:04:45
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "FWCore/Utilities/interface/Exception.h"
0012 #include "DataFormats/Provenance/interface/EventID.h"
0013 #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
0014 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0015 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
0016
0017 #include "SimGeneral/PileupInformation/interface/PileupInformation.h"
0018
0019 PileupInformation::PileupInformation(const edm::ParameterSet& config) {
0020
0021
0022 pTcut_1_ = 0.1;
0023 pTcut_2_ = 0.5;
0024 isPreMixed_ = config.getParameter<bool>("isPreMixed");
0025
0026 if (!isPreMixed_) {
0027 distanceCut_ = config.getParameter<double>("vertexDistanceCut");
0028 volumeRadius_ = config.getParameter<double>("volumeRadius");
0029 volumeZ_ = config.getParameter<double>("volumeZ");
0030 pTcut_1_ = config.getParameter<double>("pTcut_1");
0031 pTcut_2_ = config.getParameter<double>("pTcut_2");
0032
0033 PileupInfoLabel_ = consumes<PileupMixingContent>(config.getParameter<edm::InputTag>("PileupMixingLabel"));
0034
0035 PileupVtxLabel_ = consumes<PileupVertexContent>(config.getParameter<edm::InputTag>("PileupMixingLabel"));
0036
0037 LookAtTrackingTruth_ = config.getUntrackedParameter<bool>("doTrackTruth");
0038
0039 trackingTruthT_ =
0040 mayConsume<TrackingParticleCollection>(config.getParameter<edm::InputTag>("TrackingParticlesLabel"));
0041 trackingTruthV_ =
0042 mayConsume<TrackingVertexCollection>(config.getParameter<edm::InputTag>("TrackingParticlesLabel"));
0043
0044 saveVtxTimes_ = config.getParameter<bool>("saveVtxTimes");
0045
0046 MessageCategory_ = "PileupInformation";
0047
0048 edm::LogInfo(MessageCategory_) << "Setting up PileupInformation";
0049 edm::LogInfo(MessageCategory_) << "Vertex distance cut set to " << distanceCut_ << " mm";
0050 edm::LogInfo(MessageCategory_) << "Volume radius set to " << volumeRadius_ << " mm";
0051 edm::LogInfo(MessageCategory_) << "Volume Z set to " << volumeZ_ << " mm";
0052 edm::LogInfo(MessageCategory_) << "Lower pT Threshold set to " << pTcut_1_ << " GeV";
0053 edm::LogInfo(MessageCategory_) << "Upper pT Threshold set to " << pTcut_2_ << " GeV";
0054 } else {
0055 pileupSummaryToken_ =
0056 consumes<std::vector<PileupSummaryInfo> >(config.getParameter<edm::InputTag>("PileupSummaryInfoInputTag"));
0057 bunchSpacingToken_ = consumes<int>(config.getParameter<edm::InputTag>("BunchSpacingInputTag"));
0058 }
0059
0060 produces<std::vector<PileupSummaryInfo> >();
0061 produces<int>("bunchSpacing");
0062
0063
0064 }
0065
0066 void PileupInformation::produce(edm::Event& event, const edm::EventSetup& setup) {
0067 std::unique_ptr<std::vector<PileupSummaryInfo> > PSIVector(new std::vector<PileupSummaryInfo>);
0068
0069 if (isPreMixed_) {
0070 edm::Handle<std::vector<PileupSummaryInfo> > psiInput;
0071 event.getByToken(pileupSummaryToken_, psiInput);
0072
0073 std::vector<PileupSummaryInfo>::const_iterator PSiter;
0074
0075 for (PSiter = psiInput.product()->begin(); PSiter != psiInput.product()->end(); PSiter++) {
0076 PSIVector->push_back(*PSiter);
0077 }
0078
0079 edm::Handle<int> bsInput;
0080 event.getByToken(bunchSpacingToken_, bsInput);
0081 int bunchSpacing = *(bsInput.product());
0082
0083 event.put(std::move(PSIVector));
0084
0085
0086 std::unique_ptr<int> bunchSpacingP(new int(bunchSpacing));
0087 event.put(std::move(bunchSpacingP), "bunchSpacing");
0088
0089 return;
0090 }
0091
0092 edm::Handle<PileupMixingContent> MixingPileup;
0093 event.getByToken(PileupInfoLabel_, MixingPileup);
0094
0095 std::vector<int> BunchCrossings;
0096 std::vector<int> Interactions_Xing;
0097 std::vector<float> TrueInteractions_Xing;
0098 std::vector<std::vector<edm::EventID> > eventInfoList_Xing;
0099
0100 int bunchSpacing;
0101
0102
0103 const PileupMixingContent* MixInfo = MixingPileup.product();
0104 const std::vector<int>& bunchCrossing = MixInfo->getMix_bunchCrossing();
0105 const std::vector<int>& interactions = MixInfo->getMix_Ninteractions();
0106 const std::vector<float>& TrueInteractions = MixInfo->getMix_TrueInteractions();
0107 const std::vector<edm::EventID> eventInfoList = MixInfo->getMix_eventInfo();
0108
0109 bunchSpacing = MixInfo->getMix_bunchSpacing();
0110 unsigned int totalIntPU = 0;
0111
0112 for (int ib = 0; ib < (int)bunchCrossing.size(); ++ib) {
0113
0114 BunchCrossings.push_back(bunchCrossing[ib]);
0115 Interactions_Xing.push_back(interactions[ib]);
0116 TrueInteractions_Xing.push_back(TrueInteractions[ib]);
0117
0118 std::vector<edm::EventID> eventInfos;
0119 eventInfos.reserve(interactions[ib]);
0120 for (int pu = 0; pu < interactions[ib]; pu++) {
0121 eventInfos.push_back(eventInfoList[totalIntPU + pu]);
0122 }
0123 totalIntPU += (interactions[ib]);
0124 eventInfoList_Xing.push_back(eventInfos);
0125 }
0126
0127
0128 std::vector<std::vector<float> > ptHatList_Xing;
0129 std::vector<std::vector<float> > zPosList_Xing;
0130 std::vector<std::vector<float> > tPosList_Xing;
0131
0132 bool Have_pThats = false;
0133
0134 edm::Handle<PileupVertexContent> MixingPileupVtx;
0135 if (event.getByToken(PileupVtxLabel_, MixingPileupVtx)) {
0136 const PileupVertexContent* MixVtxInfo = MixingPileupVtx.product();
0137
0138 Have_pThats = true;
0139
0140 const std::vector<int>& bunchCrossing = MixInfo->getMix_bunchCrossing();
0141 const std::vector<int>& interactions = MixInfo->getMix_Ninteractions();
0142
0143 const std::vector<float>& PtHatInput = MixVtxInfo->getMix_pT_hats();
0144 const std::vector<float>& ZposInput = MixVtxInfo->getMix_z_Vtxs();
0145 const std::vector<float>& TposInput = MixVtxInfo->getMix_t_Vtxs();
0146
0147
0148
0149 unsigned int totalIntPU = 0;
0150
0151 for (int ib = 0; ib < (int)bunchCrossing.size(); ++ib) {
0152
0153
0154 std::vector<float> zposBX, tposBX;
0155 std::vector<float> pthatBX;
0156 zposBX.reserve(interactions[ib]);
0157 if (saveVtxTimes_)
0158 tposBX.reserve(interactions[ib]);
0159 pthatBX.reserve(interactions[ib]);
0160 for (int pu = 0; pu < interactions[ib]; pu++) {
0161 zposBX.push_back(ZposInput[totalIntPU + pu]);
0162 if (saveVtxTimes_)
0163 tposBX.push_back(TposInput[totalIntPU + pu]);
0164 pthatBX.push_back(PtHatInput[totalIntPU + pu]);
0165 }
0166 totalIntPU += (interactions[ib]);
0167 zPosList_Xing.push_back(zposBX);
0168 tPosList_Xing.push_back(tposBX);
0169 ptHatList_Xing.push_back(pthatBX);
0170 }
0171 }
0172
0173
0174
0175 zpositions.clear();
0176 sumpT_lowpT.clear();
0177 sumpT_highpT.clear();
0178 ntrks_lowpT.clear();
0179 ntrks_highpT.clear();
0180 event_index_.clear();
0181
0182 int lastEvent = 0;
0183
0184
0185
0186 bool HaveTrackingParticles = false;
0187
0188 edm::Handle<TrackingParticleCollection> mergedPH;
0189 edm::Handle<TrackingVertexCollection> mergedVH;
0190
0191 TrackingVertexCollection::const_iterator iVtx;
0192 TrackingVertexCollection::const_iterator iVtxTest;
0193 TrackingParticleCollection::const_iterator iTrackTest;
0194
0195 if (LookAtTrackingTruth_) {
0196 if (event.getByToken(trackingTruthT_, mergedPH) && event.getByToken(trackingTruthV_, mergedVH)) {
0197 HaveTrackingParticles = true;
0198
0199 iVtxTest = mergedVH->begin();
0200 iTrackTest = mergedPH->begin();
0201 }
0202 }
0203
0204 int nminb_vtx = 0;
0205
0206
0207
0208 std::vector<int>::iterator BXIter;
0209 std::vector<int>::iterator InteractionsIter = Interactions_Xing.begin();
0210 std::vector<float>::iterator TInteractionsIter = TrueInteractions_Xing.begin();
0211 std::vector<std::vector<edm::EventID> >::iterator TEventInfoIter = eventInfoList_Xing.begin();
0212
0213 std::vector<std::vector<float> >::iterator zPosIter;
0214 std::vector<std::vector<float> >::iterator tPosIter;
0215 std::vector<std::vector<float> >::iterator pThatIter;
0216
0217 if (Have_pThats) {
0218 zPosIter = zPosList_Xing.begin();
0219 tPosIter = tPosList_Xing.begin();
0220 pThatIter = ptHatList_Xing.begin();
0221 }
0222
0223
0224
0225 for (BXIter = BunchCrossings.begin(); BXIter != BunchCrossings.end();
0226 ++BXIter, ++InteractionsIter, ++TInteractionsIter, ++TEventInfoIter) {
0227
0228
0229 if (HaveTrackingParticles) {
0230
0231 for (iVtx = iVtxTest; iVtx != mergedVH->end(); ++iVtx) {
0232 if (iVtx->eventId().bunchCrossing() == (*BXIter)) {
0233
0234 if (iVtx->eventId().event() != lastEvent) {
0235
0236
0237 float zpos = 0.;
0238 zpos = iVtx->position().z();
0239 zpositions.push_back(zpos);
0240 sumpT_lowpT.push_back(0.);
0241 sumpT_highpT.push_back(0.);
0242 ntrks_lowpT.push_back(0);
0243 ntrks_highpT.push_back(0);
0244
0245 lastEvent = iVtx->eventId().event();
0246 iVtxTest = --iVtx;
0247
0248
0249
0250 event_index_.insert(myindex::value_type(lastEvent, nminb_vtx));
0251
0252 ++nminb_vtx;
0253
0254 continue;
0255 }
0256 }
0257 }
0258
0259
0260
0261 for (TrackingParticleCollection::const_iterator iTrack = iTrackTest; iTrack != mergedPH->end(); ++iTrack) {
0262 bool FoundTrk = false;
0263
0264 float zpos = 0.;
0265
0266 if (iTrack->eventId().bunchCrossing() == (*BXIter) && iTrack->eventId().event() > 0) {
0267 FoundTrk = true;
0268 int correct_index = event_index_[iTrack->eventId().event()];
0269
0270
0271
0272 zpos = zpositions[correct_index];
0273 if (iTrack->matchedHit() > 0) {
0274 if (fabs(iTrack->parentVertex()->position().z() - zpos) <
0275 0.1) {
0276
0277 float Tpx = iTrack->p4().px();
0278 float Tpy = iTrack->p4().py();
0279 float TpT = sqrt(Tpx * Tpx + Tpy * Tpy);
0280 if (TpT > pTcut_1_) {
0281 sumpT_lowpT[correct_index] += TpT;
0282 ++ntrks_lowpT[correct_index];
0283 }
0284 if (TpT > pTcut_2_) {
0285 sumpT_highpT[correct_index] += TpT;
0286 ++ntrks_highpT[correct_index];
0287 }
0288 }
0289 }
0290 } else {
0291 if (FoundTrk) {
0292 iTrackTest = --iTrack;
0293 --iTrackTest;
0294 break;
0295 }
0296 }
0297
0298 }
0299
0300 }
0301
0302
0303
0304
0305
0306
0307 if (!HaveTrackingParticles) {
0308
0309 zpositions.push_back(-999.);
0310 sumpT_lowpT.push_back(0.);
0311 sumpT_highpT.push_back(0.);
0312 ntrks_lowpT.push_back(0);
0313 ntrks_highpT.push_back(0);
0314 }
0315
0316 if (Have_pThats) {
0317 PileupSummaryInfo PSI_bunch = PileupSummaryInfo((*InteractionsIter),
0318 (*zPosIter),
0319 (*tPosIter),
0320 sumpT_lowpT,
0321 sumpT_highpT,
0322 ntrks_lowpT,
0323 ntrks_highpT,
0324 (*TEventInfoIter),
0325 (*pThatIter),
0326 (*BXIter),
0327 (*TInteractionsIter),
0328 bunchSpacing);
0329 PSIVector->push_back(PSI_bunch);
0330
0331 zPosIter++;
0332 tPosIter++;
0333 pThatIter++;
0334 } else {
0335 std::vector<float> zposZeros((*TEventInfoIter).size(), 0);
0336 std::vector<float> tposZeros((*TEventInfoIter).size(), 0);
0337 std::vector<float> pThatZeros((*TEventInfoIter).size(), 0);
0338
0339 PileupSummaryInfo PSI_bunch = PileupSummaryInfo((*InteractionsIter),
0340 zposZeros,
0341 tposZeros,
0342 sumpT_lowpT,
0343 sumpT_highpT,
0344 ntrks_lowpT,
0345 ntrks_highpT,
0346 (*TEventInfoIter),
0347 pThatZeros,
0348 (*BXIter),
0349 (*TInteractionsIter),
0350 bunchSpacing);
0351
0352 PSIVector->push_back(PSI_bunch);
0353 }
0354
0355
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365
0366
0367
0368
0369 event_index_.clear();
0370 zpositions.clear();
0371 sumpT_lowpT.clear();
0372 sumpT_highpT.clear();
0373 ntrks_lowpT.clear();
0374 ntrks_highpT.clear();
0375 nminb_vtx = 0;
0376 lastEvent = 0;
0377
0378 }
0379
0380
0381
0382 event.put(std::move(PSIVector));
0383
0384
0385 std::unique_ptr<int> bunchSpacingP(new int(bunchSpacing));
0386 event.put(std::move(bunchSpacingP), "bunchSpacing");
0387 }
0388
0389 DEFINE_FWK_MODULE(PileupInformation);