Go to the documentation of this file.00001
00002
00003
00004
00005
00006
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
00023
00024 pTcut_1_ = 0.1;
00025 pTcut_2_ = 0.5;
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
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;
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) {
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
00071 BunchCrossings.push_back(bunchCrossing[ib]);
00072 Interactions_Xing.push_back(interactions[ib]);
00073 }
00074 }
00075 else{
00076
00077
00078 edm::Handle<CrossingFrame<SimVertex> > cfSimVertexes;
00079 event.getByLabel("mix", simHitLabel_, cfSimVertexes);
00080
00081
00082 simVertexes_ = std::auto_ptr<MixCollection<SimVertex> >( new MixCollection<SimVertex>(cfSimVertexes.product()) );
00083
00084 int index = 0;
00085
00086
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
00098 for (MixCollection<SimVertex>::MixItr iterator = simVertexes_->begin(); iterator != simVertexes_->end(); ++iterator, ++index)
00099 {
00100
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
00130 BunchCrossings.push_back((*viter));
00131 Interactions_Xing.push_back((*liter));
00132 }
00133 }
00134
00135
00136
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;
00153
00154 int lastBunchCrossing = 0;
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
00162
00163
00164 std::vector<int>::iterator BXIter;
00165 std::vector<int>::iterator InteractionsIter = Interactions_Xing.begin();
00166
00167
00168
00169 for( BXIter = BunchCrossings.begin(); BXIter != BunchCrossings.end(); ++BXIter, ++InteractionsIter) {
00170
00171
00172
00173 for (iVtx = iVtxTest; iVtx != mergedVH->end(); ++iVtx) {
00174
00175 if(iVtx->eventId().bunchCrossing() == (*BXIter) ) {
00176
00177 if(iVtx->eventId().event() != lastEvent) {
00178
00179
00180
00181 float zpos = 0.;
00182 zpos = iVtx->position().z();
00183 zpositions.push_back(zpos);
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;
00191
00192
00193
00194 event_index_.insert(myindex::value_type(lastEvent,nminb_vtx));
00195
00196 ++nminb_vtx;
00197
00198 continue;
00199 }
00200 }
00201 }
00202
00203
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
00217
00218 zpos = zpositions[correct_index];
00219 if(iTrack->matchedHit()>0) {
00220 if(fabs(iTrack->parentVertex()->position().z()-zpos)<0.1) {
00221
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;
00240 --iTrackTest;
00241 break;
00242 }
00243
00244 }
00245
00246 }
00247
00248
00249
00250
00251
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
00264
00265
00266
00267
00268
00269
00270
00271
00272
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 }
00290
00291
00292
00293 event.put(PSIVector);
00294
00295
00296 }
00297
00298
00299 DEFINE_FWK_MODULE(PileupInformation);