CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/SimGeneral/PileupInformation/plugins/PileupInformation.cc

Go to the documentation of this file.
00001 // File: PileupInformation.cc
00002 // Description:  adds pileup information object to event
00003 // Author:  Mike Hildreth
00004 //
00005 // Adds a vector of PileupSummaryInfo objects to the event. 
00006 // One for each bunch crossing.
00007 //
00008 //--------------------------------------------
00009 
00010 #include "FWCore/Framework/interface/MakerMacros.h"
00011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00012 
00013 #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
00014 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
00015 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
00016 
00017 #include "SimGeneral/PileupInformation/interface/PileupInformation.h"
00018 #include "SimDataFormats/PileupSummaryInfo/interface/PileupMixingContent.h"
00019 
00020 PileupInformation::PileupInformation(const edm::ParameterSet & config) 
00021 {
00022     // Initialize global parameters
00023 
00024     pTcut_1_                = 0.1;
00025     pTcut_2_                = 0.5; // defaults                                                       
00026     distanceCut_            = config.getParameter<double>("vertexDistanceCut");
00027     volumeRadius_           = config.getParameter<double>("volumeRadius");
00028     volumeZ_                = config.getParameter<double>("volumeZ");
00029     pTcut_1_                = config.getParameter<double>("pTcut_1");
00030     pTcut_2_                = config.getParameter<double>("pTcut_2");
00031 
00032     trackingTruth_  = config.getParameter<std::string>("TrackingParticlesLabel");
00033 
00034     simHitLabel_            = config.getParameter<std::string>("simHitLabel");
00035 
00036     MessageCategory_       = "PileupInformation";
00037 
00038     edm::LogInfo (MessageCategory_) << "Setting up PileupInformation";
00039     edm::LogInfo (MessageCategory_) << "Vertex distance cut set to " << distanceCut_  << " mm";
00040     edm::LogInfo (MessageCategory_) << "Volume radius set to "       << volumeRadius_ << " mm";
00041     edm::LogInfo (MessageCategory_) << "Volume Z      set to "       << volumeZ_      << " mm";
00042     edm::LogInfo (MessageCategory_) << "Lower pT Threshold set to "       << pTcut_1_      << " GeV";
00043     edm::LogInfo (MessageCategory_) << "Upper pT Threshold set to "       << pTcut_2_      << " GeV";
00044 
00045 
00046     produces< std::vector<PileupSummaryInfo> >();
00047     //produces<PileupSummaryInfo>();
00048 }
00049 
00050 
00051 void PileupInformation::produce(edm::Event &event, const edm::EventSetup & setup)
00052 {
00053 
00054   std::auto_ptr<std::vector<PileupSummaryInfo> > PSIVector(new std::vector<PileupSummaryInfo>);
00055 
00056   edm::Handle< PileupMixingContent > MixingPileup;  // Get True pileup information from MixingModule
00057   event.getByLabel("mix", MixingPileup);
00058 
00059   std::vector<int> BunchCrossings;
00060   std::vector<int> Interactions_Xing;
00061 
00062   const PileupMixingContent* MixInfo = MixingPileup.product();
00063 
00064   if(MixInfo) {  // extract information - way easier than counting vertices
00065 
00066     const std::vector<int> bunchCrossing = MixInfo->getMix_bunchCrossing();
00067     const std::vector<int> interactions = MixInfo->getMix_Ninteractions();
00068 
00069     for(int ib=0; ib<(int)bunchCrossing.size(); ++ib){
00070       //      std::cout << " bcr, nint " << bunchCrossing[ib] << " " << interactions[ib] << std::endl;
00071       BunchCrossings.push_back(bunchCrossing[ib]);
00072       Interactions_Xing.push_back(interactions[ib]);
00073     }
00074   }
00075   else{ //If no Mixing Truth, work with SimVertices  (probably should throw an exception, but...)
00076 
00077     // Collect all the simvertex from the crossing frame                                            
00078     edm::Handle<CrossingFrame<SimVertex> > cfSimVertexes;
00079     event.getByLabel("mix", simHitLabel_, cfSimVertexes);
00080 
00081     // Create a mix collection from one simvertex collection                                        
00082     simVertexes_ = std::auto_ptr<MixCollection<SimVertex> >( new MixCollection<SimVertex>(cfSimVertexes.product()) );
00083 
00084     int index = 0;
00085     // Solution to the problem of not having vertexId
00086     //    bool FirstL = true;
00087     EncodedEventIdToIndex vertexId;
00088     EncodedEventId oldEventId;
00089     unsigned int oldVertexId = 0;
00090     int oldBX = -1000;
00091     int oldEvent = 0;
00092 
00093     std::vector<int> BunchCrossings2;
00094     std::list<int> Interactions_Xing2;
00095 
00096 
00097     // Loop for finding repeated vertexId (vertexId problem hack)                                   
00098     for (MixCollection<SimVertex>::MixItr iterator = simVertexes_->begin(); iterator != simVertexes_->end(); ++iterator, ++index)
00099       {
00100         //      std::cout << " SimVtx eventid, vertexid " << iterator->eventId().event() << " " << iterator->eventId().bunchCrossing() << std::endl;
00101         if (!index || iterator->eventId() != oldEventId)
00102           {
00103             if(iterator->eventId().bunchCrossing()==0 && iterator->eventId().event()==0){
00104               continue;
00105             }
00106             if(iterator->eventId().bunchCrossing() != oldBX) {
00107               BunchCrossings2.push_back(iterator->eventId().bunchCrossing());
00108               Interactions_Xing2.push_back(iterator->eventId().event());
00109               oldBX = iterator->eventId().bunchCrossing();
00110               oldEvent = iterator->eventId().event();
00111             }
00112             else { Interactions_Xing2.pop_back();
00113               Interactions_Xing2.push_back(iterator->eventId().event());
00114               oldEvent = iterator->eventId().event();
00115             }
00116 
00117 
00118             oldEventId = iterator->eventId();
00119             oldVertexId = iterator->vertexId();
00120             continue;
00121           }
00122 
00123       }
00124 
00125     std::vector<int>::iterator viter;
00126     std::list<int>::iterator liter = Interactions_Xing2.begin();
00127 
00128     for(viter = BunchCrossings2.begin(); viter != BunchCrossings2.end(); ++viter, ++liter){
00129       //std::cout << " bcr, nint from VTX " << (*viter) << " " << (*liter) << std::endl;
00130       BunchCrossings.push_back((*viter));
00131       Interactions_Xing.push_back((*liter));
00132     }
00133   } // end of look at SimVertices
00134 
00135 
00136   //Now, get information on valid particles that look like they could be in the tracking volume
00137 
00138   edm::Handle<TrackingParticleCollection> mergedPH;
00139   edm::Handle<TrackingVertexCollection>   mergedVH;
00140 
00141   event.getByLabel(trackingTruth_, mergedPH);
00142   event.getByLabel(trackingTruth_, mergedVH);
00143 
00144 
00145   zpositions.clear();
00146   sumpT_lowpT.clear();
00147   sumpT_highpT.clear();
00148   ntrks_lowpT.clear();
00149   ntrks_highpT.clear();
00150   event_index_.clear();
00151 
00152   int lastEvent = 0; // zero is the true MC hard-scatter event
00153 
00154   int lastBunchCrossing = 0; // 0 is the true bunch crossing, should always come first.
00155 
00156   TrackingVertexCollection::const_iterator iVtx;
00157   TrackingVertexCollection::const_iterator iVtxTest = mergedVH->begin();
00158   TrackingParticleCollection::const_iterator iTrackTest = mergedPH->begin();
00159 
00160   int nminb_vtx = 0;
00161   //  bool First = true;
00162   //  bool flag_new = false;
00163 
00164   std::vector<int>::iterator BXIter;
00165   std::vector<int>::iterator InteractionsIter = Interactions_Xing.begin();
00166 
00167   // loop over the bunch crossings and interactions we have extracted 
00168 
00169   for( BXIter = BunchCrossings.begin(); BXIter != BunchCrossings.end(); ++BXIter, ++InteractionsIter) {
00170 
00171     //std::cout << "looking for BX: " << (*BXIter) << std::endl;
00172 
00173     for (iVtx = iVtxTest; iVtx != mergedVH->end(); ++iVtx) {     
00174 
00175       if(iVtx->eventId().bunchCrossing() == (*BXIter) ) { // found first vertex in this bunch crossing
00176 
00177         if(iVtx->eventId().event() != lastEvent) {
00178 
00179           //std::cout << "BX,event " << iVtx->eventId().bunchCrossing() << " " << iVtx->eventId().event() << std::endl;
00180 
00181              float zpos = 0.;
00182              zpos = iVtx->position().z();
00183              zpositions.push_back(zpos);  //save z position of each vertex
00184              sumpT_lowpT.push_back(0.);
00185              sumpT_highpT.push_back(0.);
00186              ntrks_lowpT.push_back(0);
00187              ntrks_highpT.push_back(0);
00188 
00189              lastEvent = iVtx->eventId().event();
00190              iVtxTest = --iVtx; // just for security
00191 
00192              // turns out events aren't sequential... save map of indices
00193 
00194              event_index_.insert(myindex::value_type(lastEvent,nminb_vtx));
00195              
00196              ++nminb_vtx;
00197 
00198              continue;
00199         }
00200       }
00201     }
00202     
00203     // next loop over tracks to get information
00204 
00205     for (TrackingParticleCollection::const_iterator iTrack = iTrackTest; iTrack != mergedPH->end(); ++iTrack)
00206       {
00207         bool FoundTrk = false;
00208 
00209         float zpos=0.;
00210 
00211         if(iTrack->eventId().bunchCrossing() == (*BXIter) && iTrack->eventId().event() > 0 )
00212           {
00213             FoundTrk = true;
00214             int correct_index = event_index_[iTrack->eventId().event()];
00215 
00216             //std::cout << " track index, correct index " << iTrack->eventId().event() << " " << correct_index << std::endl;
00217 
00218             zpos = zpositions[correct_index];
00219             if(iTrack->matchedHit()>0) {
00220               if(fabs(iTrack->parentVertex()->position().z()-zpos)<0.1) {  //make sure track really comes from this vertex
00221                 //std::cout << *iTrack << std::endl;                                              
00222                 float Tpx = iTrack->p4().px();
00223                 float Tpy = iTrack->p4().py();
00224                 float TpT = sqrt(Tpx*Tpx + Tpy*Tpy);
00225                 if( TpT>pTcut_1_ ) {
00226                   sumpT_lowpT[correct_index]+=TpT;
00227                   ++ntrks_lowpT[correct_index];
00228                 }
00229                 if( TpT>pTcut_2_ ){
00230                   sumpT_highpT[correct_index]+=TpT;
00231                   ++ntrks_highpT[correct_index];
00232                 }
00233               }
00234             }
00235           }
00236         else{
00237           if(FoundTrk) {
00238 
00239             iTrackTest = --iTrack;  // reset so we can start over next time
00240             --iTrackTest;  // just to be sure
00241             break;
00242           }
00243         
00244         }
00245         
00246       } // end of track loop
00247 
00248         // now that we have all of the track information for a given bunch crossing, 
00249         // make PileupSummary for this one and move on
00250 
00251         //        std::cout << "Making PSI for bunch " << lastBunchCrossing << std::endl;
00252 
00253         PileupSummaryInfo       PSI_bunch = PileupSummaryInfo(
00254                                                               (*InteractionsIter),
00255                                                               zpositions,
00256                                                               sumpT_lowpT,
00257                                                               sumpT_highpT,
00258                                                               ntrks_lowpT,
00259                                                               ntrks_highpT,
00260                                                               (*BXIter)
00261                                                               );
00262 
00263         //std::cout << " " << std::endl;
00264         //std::cout << "Adding Bunch Crossing, nint " << (*BXIter) << " " <<  (*InteractionsIter) << std::endl;
00265  
00266         // for(int iv = 0; iv<nminb_vtx; ++iv){
00267             
00268           //  std::cout << "Z position " << zpositions[iv] << std::endl;
00269           //  std::cout << "ntrks_lowpT " << ntrks_lowpT[iv] << std::endl;
00270           //  std::cout << "sumpT_lowpT " << sumpT_lowpT[iv] << std::endl;
00271           //  std::cout << "ntrks_highpT " << ntrks_highpT[iv] << std::endl;
00272           //  std::cout << "sumpT_highpT " << sumpT_highpT[iv] << std::endl;
00273           // }
00274 
00275         PSIVector->push_back(PSI_bunch);
00276 
00277         lastBunchCrossing = iVtx->eventId().bunchCrossing();
00278 
00279         event_index_.clear();
00280         zpositions.clear();
00281         sumpT_lowpT.clear();
00282         sumpT_highpT.clear();
00283         ntrks_lowpT.clear();
00284         ntrks_highpT.clear();
00285         nminb_vtx = 0;
00286         lastEvent=0;
00287 
00288 
00289   } // end of loop over bunch crossings
00290 
00291   // put our vector of PileupSummaryInfo objects into the event.
00292 
00293   event.put(PSIVector);
00294 
00295 
00296 }
00297 
00298 
00299 DEFINE_FWK_MODULE(PileupInformation);