Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 10:04:45

0001 // File: PileupInformation.cc
0002 // Description:  adds pileup information object to event
0003 // Author:  Mike Hildreth
0004 //
0005 // Adds a vector of PileupSummaryInfo objects to the event.
0006 // One for each bunch crossing.
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   // Initialize global parameters
0021 
0022   pTcut_1_ = 0.1;
0023   pTcut_2_ = 0.5;  // defaults
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   //produces<PileupSummaryInfo>();
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     //add bunch spacing to the event as a seperate integer for use by downstream modules
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;  // Get True pileup information from MixingModule
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   // extract information - way easier than counting vertices
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     //      std::cout << " bcr, nint " << bunchCrossing[ib] << " " << interactions[ib] << std::endl;
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   // store information from pileup vertices, if it's in the event. Have to loop on interactions again.
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;          // Get True pileup information from MixingModule
0135   if (event.getByToken(PileupVtxLabel_, MixingPileupVtx)) {  // extract information - way easier than counting vertices
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     // store information from pileup vertices, if it's in the event:
0148 
0149     unsigned int totalIntPU = 0;
0150 
0151     for (int ib = 0; ib < (int)bunchCrossing.size(); ++ib) {
0152       //      std::cout << " bcr, nint " << bunchCrossing[ib] << " " << interactions[ib] << std::endl;
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   }  // end of VertexInfo block
0172 
0173   //Now, get information on valid particles that look like they could be in the tracking volume
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;  // zero is the true MC hard-scatter event
0183 
0184   // int lastBunchCrossing = 0; // 0 is the true bunch crossing, should always come first.
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   //  bool First = true;
0206   //  bool flag_new = false;
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   // loop over the bunch crossings and interactions we have extracted
0224 
0225   for (BXIter = BunchCrossings.begin(); BXIter != BunchCrossings.end();
0226        ++BXIter, ++InteractionsIter, ++TInteractionsIter, ++TEventInfoIter) {
0227     //std::cout << "looking for BX: " << (*BXIter) << std::endl;
0228 
0229     if (HaveTrackingParticles) {  // leave open the option of not producing TrackingParticles and just keeping # interactions
0230 
0231       for (iVtx = iVtxTest; iVtx != mergedVH->end(); ++iVtx) {
0232         if (iVtx->eventId().bunchCrossing() == (*BXIter)) {  // found first vertex in this bunch crossing
0233 
0234           if (iVtx->eventId().event() != lastEvent) {
0235             //std::cout << "BX,event " << iVtx->eventId().bunchCrossing() << " " << iVtx->eventId().event() << std::endl;
0236 
0237             float zpos = 0.;
0238             zpos = iVtx->position().z();
0239             zpositions.push_back(zpos);  //save z position of each vertex
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;  // just for security
0247 
0248             // turns out events aren't sequential... save map of indices
0249 
0250             event_index_.insert(myindex::value_type(lastEvent, nminb_vtx));
0251 
0252             ++nminb_vtx;
0253 
0254             continue;
0255           }
0256         }
0257       }
0258 
0259       // next loop over tracks to get information
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           //std::cout << " track index, correct index " << iTrack->eventId().event() << " " << correct_index << std::endl;
0271 
0272           zpos = zpositions[correct_index];
0273           if (iTrack->matchedHit() > 0) {
0274             if (fabs(iTrack->parentVertex()->position().z() - zpos) <
0275                 0.1) {  //make sure track really comes from this vertex
0276               //std::cout << *iTrack << std::endl;
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;  // reset so we can start over next time
0293             --iTrackTest;           // just to be sure
0294             break;
0295           }
0296         }
0297 
0298       }  // end of track loop
0299 
0300     }  // end of check that we have TrackingParticles to begin with...
0301 
0302     // now that we have all of the track information for a given bunch crossing,
0303     // make PileupSummary for this one and move on
0304 
0305     //    std::cout << "Making PSI for bunch " << lastBunchCrossing << std::endl;
0306 
0307     if (!HaveTrackingParticles) {  // stick in one value so we don't have empty vectors
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     //std::cout << " " << std::endl;
0355     //std::cout << "Adding Bunch Crossing, nint " << (*BXIter) << " " <<  (*InteractionsIter) << std::endl;
0356 
0357     //for(int iv = 0; iv<(*InteractionsIter); ++iv){
0358 
0359     // std::cout << "Z position " << zpositions[iv] << std::endl;
0360     // std::cout << "ntrks_lowpT " << ntrks_lowpT[iv] << std::endl;
0361     // std::cout << "sumpT_lowpT " << sumpT_lowpT[iv] << std::endl;
0362     // std::cout << "ntrks_highpT " << ntrks_highpT[iv] << std::endl;
0363     // std::cout << "sumpT_highpT " << sumpT_highpT[iv] << std::endl;
0364     //std::cout << iv << " " << PSI_bunch.getPU_EventID()[iv] << std::endl;
0365     //}
0366 
0367     // if(HaveTrackingParticles) lastBunchCrossing = iVtx->eventId().bunchCrossing();
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   }  // end of loop over bunch crossings
0379 
0380   // put our vector of PileupSummaryInfo objects into the event.
0381 
0382   event.put(std::move(PSIVector));
0383 
0384   //add bunch spacing to the event as a seperate integer for use by downstream modules
0385   std::unique_ptr<int> bunchSpacingP(new int(bunchSpacing));
0386   event.put(std::move(bunchSpacingP), "bunchSpacing");
0387 }
0388 
0389 DEFINE_FWK_MODULE(PileupInformation);