Go to the documentation of this file.00001
00002
00003
00004 #include "PhysicsTools/HepMCCandAlgos/interface/FlavorHistoryProducer.h"
00005
00006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00007 #include "DataFormats/Math/interface/deltaR.h"
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 using namespace std;
00035 using namespace reco;
00036 using namespace edm;
00037
00038 FlavorHistoryProducer::FlavorHistoryProducer( const ParameterSet & p ) :
00039 src_( p.getParameter<InputTag>( "src" ) ),
00040 matchedSrc_( p.getParameter<InputTag>( "matchedSrc") ),
00041 matchDR_ ( p.getParameter<double> ("matchDR") ),
00042 pdgIdToSelect_( p.getParameter<int> ("pdgIdToSelect") ),
00043 ptMinParticle_( p.getParameter<double>( "ptMinParticle") ),
00044 ptMinShower_( p.getParameter<double>( "ptMinShower") ),
00045 etaMaxParticle_( p.getParameter<double>( "etaMaxParticle" )),
00046 etaMaxShower_( p.getParameter<double>( "etaMaxShower" )),
00047 flavorHistoryName_( p.getParameter<string>("flavorHistoryName") ),
00048 verbose_( p.getUntrackedParameter<bool>( "verbose" ) )
00049 {
00050 produces<FlavorHistoryEvent >(flavorHistoryName_);
00051 }
00052
00053 FlavorHistoryProducer::~FlavorHistoryProducer() {
00054 }
00055
00056 void FlavorHistoryProducer::produce( Event& evt, const EventSetup& )
00057 {
00058 if ( verbose_ ) cout << "---------- Hello from FlavorHistoryProducer! -----" << endl;
00059
00060
00061 Handle<CandidateView > particlesViewH;
00062 evt.getByLabel( src_, particlesViewH );
00063
00064 Handle<CandidateView> matchedView;
00065 evt.getByLabel( matchedSrc_, matchedView );
00066
00067
00068 vector<const Candidate* > particles;
00069 for( CandidateView::const_iterator p = particlesViewH->begin(); p != particlesViewH->end(); ++p ) {
00070 particles.push_back(&*p);
00071 }
00072
00073
00074 vector<int> partonIndices;
00075
00076 vector<int> progenitorIndices;
00077
00078 vector<int> sisterIndices;
00079
00080 vector<FlavorHistory::FLAVOR_T> flavorSources;
00081
00082
00083 auto_ptr<FlavorHistoryEvent > flavorHistoryEvent ( new FlavorHistoryEvent () ) ;
00084
00085
00086
00087
00088 vector<const Candidate* >::size_type j;
00089 vector<const Candidate* >::size_type j_max = particles.size();
00090 for( j=0; j<j_max; ++j ) {
00091
00092
00093 const Candidate *p = particles[j];
00094
00095 vector<Candidate const *>::size_type partonIndex=j;
00096 vector<Candidate const *>::size_type progenitorIndex=j_max;
00097 bool foundProgenitor = false;
00098 FlavorHistory::FLAVOR_T flavorSource=FlavorHistory::FLAVOR_NULL;
00099
00100
00101 int idabs = std::abs( (p)->pdgId() );
00102 int nDa = (p)->numberOfDaughters();
00103
00104
00105
00106 if ( idabs <= FlavorHistory::tQuarkId && p->status() == 2 &&
00107 idabs == pdgIdToSelect_ ) {
00108
00109 if ( nDa > 0 ) {
00110
00111 if((p)->pt() > ptMinShower_ && fabs((p)->eta())<etaMaxShower_) {
00112
00113 if(verbose_) cout << "--------------------------" << endl;
00114 if(verbose_) cout << "Processing particle " << j << ", status = " << p->status() << ", pdg id = " << p->pdgId() << endl;
00115
00116
00117
00118 vector<Candidate const *> allParents;
00119 getAncestors( *p, allParents );
00120
00121 if(verbose_) {
00122 cout << "Parents : " << endl;
00123 vector<Candidate const *>::const_iterator iprint = allParents.begin(),
00124 iprintend = allParents.end();
00125 for ( ; iprint != iprintend; ++iprint )
00126 cout << " status = " << (*iprint)->status() << ", pdg id = " << (*iprint)->pdgId() << ", pt = " << (*iprint)->pt() << endl;
00127 }
00128
00129
00130
00131 bool status3AncestorOfSameFlavor = false;
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153 vector<Candidate const *>::size_type a_size = allParents.size();
00154
00155
00156
00157
00158 for( vector<Candidate const *>::size_type i=0 ; i < a_size && !foundProgenitor; ++i ) {
00159 const Candidate * aParent=allParents[i];
00160 vector<Candidate const *>::const_iterator found = find(particles.begin(),particles.end(),aParent);
00161
00162
00163
00164 progenitorIndex = found - particles.begin();
00165
00166 int aParentId = std::abs(aParent->pdgId());
00167
00168
00169 if ( aParent->status() == 3 && aParent->pdgId() == p->pdgId() ) status3AncestorOfSameFlavor = true;
00170
00171
00172
00173
00174 if( progenitorIndex > 5 ) {
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184 if( aParent->numberOfMothers() == 2 &&
00185 aParent->pdgId() == p->pdgId() && aParent->status() == 3 ) {
00186 if(verbose_) cout << "Matrix Element progenitor" << endl;
00187 if ( flavorSource == FlavorHistory::FLAVOR_NULL ) flavorSource = FlavorHistory::FLAVOR_ME;
00188 foundProgenitor = true;
00189 }
00190
00191 else if( (aParentId>pdgIdToSelect_ && aParentId<FlavorHistory::gluonId) ||
00192 aParentId > FlavorHistory::gluonId ) {
00193 if(verbose_) cout << "Flavor decay progenitor" << endl;
00194 if ( flavorSource == FlavorHistory::FLAVOR_NULL ) flavorSource = FlavorHistory::FLAVOR_DECAY;
00195 foundProgenitor = true;
00196 }
00197 }
00198 }
00199
00200
00201
00202 if ( !foundProgenitor ) {
00203 if ( flavorSource == FlavorHistory::FLAVOR_NULL ) flavorSource = FlavorHistory::FLAVOR_GS;
00204
00205 for( vector<Candidate const *>::size_type i=0 ; i < a_size && !foundProgenitor; ++i ) {
00206 const Candidate * aParent=allParents[i];
00207 vector<Candidate const *>::const_iterator found = find(particles.begin(),particles.end(),aParent);
00208
00209 progenitorIndex = found - particles.begin();
00210
00211 if ( aParent->numberOfMothers() == 1 && aParent->status() == 3 ) {
00212 foundProgenitor = true;
00213 }
00214 }
00215 }
00216
00217 }
00218 }
00219
00220
00221
00222
00223 if ( !foundProgenitor ) progenitorIndex = j_max;
00224
00225
00226
00227 partonIndices.push_back( partonIndex );
00228 progenitorIndices.push_back( progenitorIndex );
00229 flavorSources.push_back(flavorSource);
00230 sisterIndices.push_back( -1 );
00231
00232 }
00233 }
00234
00235
00236
00237
00238
00239 if ( verbose_ ) cout << "Making sisters" << endl;
00240
00241 if ( partonIndices.size() == progenitorIndices.size() && partonIndices.size() > 0 ) {
00242
00243 for ( unsigned int ii = 0; ii < partonIndices.size(); ++ii ) {
00244
00245 const Candidate * candi = particles[partonIndices[ii]];
00246 if ( verbose_ ) cout << " Sister candidate particle 1: " << ii << ", pdgid = " << candi->pdgId() << ", status = " << candi->status() << endl;
00247
00248
00249
00250 for ( unsigned int jj = 0; jj < partonIndices.size(); ++jj ) {
00251 if ( ii != jj ) {
00252 const Candidate * candj = particles[partonIndices[jj]];
00253 if ( verbose_ ) cout << " Sister candidate particle 2: " << jj << ", pdgid = " << candj->pdgId() << ", status = " << candj->status() << endl;
00254
00255
00256 if ( candi->pdgId() == -1 * candj->pdgId() && candi->status() == candj->status() ) {
00257 sisterIndices[ii] = partonIndices[jj];
00258 if ( verbose_ ) cout << "Parton " << partonIndices[ii] << " has sister " << sisterIndices[ii] << endl;
00259 }
00260 }
00261 }
00262
00263
00264 if ( sisterIndices[ii] < 0 ) {
00265 if ( verbose_ ) cout << "No sister. Classified as flavor excitation" << endl;
00266 flavorSources[ii] = FlavorHistory::FLAVOR_EXC;
00267 }
00268
00269 if ( verbose_ ) cout << "Getting matched jet" << endl;
00270
00271 CandidateView::const_iterator matched = getClosestJet( matchedView, *candi );
00272 CandidateView::const_iterator matchedBegin = matchedView->begin();
00273 CandidateView::const_iterator matchedEnd = matchedView->end();
00274
00275
00276 if ( verbose_ ) cout << "Getting sister jet" << endl;
00277
00278 CandidateView::const_iterator sister =
00279 ( sisterIndices[ii] >= 0 && static_cast<unsigned int>(sisterIndices[ii]) < particles.size() ) ?
00280 getClosestJet( matchedView, *particles[sisterIndices[ii]] ) :
00281 matchedEnd ;
00282
00283 if ( verbose_ ) cout << "Making matched shallow clones : ";
00284 ShallowClonePtrCandidate matchedCand ;
00285 if ( matched != matchedEnd ) {
00286 if ( verbose_ ) cout << " found" << endl;
00287 matchedCand = ShallowClonePtrCandidate( CandidatePtr(matchedView, matched - matchedBegin ) );
00288 } else {
00289 if ( verbose_ ) cout << " NOT found" << endl;
00290 }
00291
00292 if ( verbose_ ) cout << "Making sister shallow clones : ";
00293 ShallowClonePtrCandidate sisterCand;
00294 if ( sister != matchedEnd ) {
00295 if ( verbose_ ) cout << " found" << endl;
00296 sisterCand = ShallowClonePtrCandidate( CandidatePtr(matchedView, sister - matchedBegin ) );
00297 } else {
00298 if ( verbose_ ) cout << " NOT found" << endl;
00299 }
00300
00301 if ( verbose_ ) cout << "Making history object" << endl;
00302
00303 FlavorHistory history (flavorSources[ii],
00304 particlesViewH,
00305 partonIndices[ii], progenitorIndices[ii], sisterIndices[ii],
00306 matchedCand,
00307 sisterCand );
00308 if ( verbose_ ) cout << "Adding flavor history : " << history << endl;
00309 flavorHistoryEvent->push_back( history );
00310 }
00311 }
00312
00313
00314 if ( flavorHistoryEvent->size() > 0 ) {
00315
00316
00317 if ( verbose_ ) cout << "Caching flavor history event" << endl;
00318 flavorHistoryEvent->cache();
00319
00320 if ( verbose_ ) {
00321 cout << "Outputting pdg id = " << pdgIdToSelect_ << " with nelements = " << flavorHistoryEvent->size() << endl;
00322 vector<FlavorHistory>::const_iterator i = flavorHistoryEvent->begin(),
00323 iend = flavorHistoryEvent->end();
00324 for ( ; i !=iend; ++i ) {
00325 cout << *i << endl;
00326 }
00327 }
00328 }
00329
00330
00331 evt.put( flavorHistoryEvent, flavorHistoryName_ );
00332 }
00333
00334
00335
00336 void FlavorHistoryProducer::getAncestors(const Candidate &c,
00337 vector<Candidate const *> & moms )
00338 {
00339
00340 if( c.numberOfMothers() == 1 ) {
00341 const Candidate * dau = &c;
00342 const Candidate * mom = c.mother();
00343 while ( dau->numberOfMothers() != 0) {
00344 moms.push_back( dau );
00345 dau = mom ;
00346 mom = dau->mother();
00347 }
00348 }
00349 }
00350
00351
00352 CandidateView::const_iterator
00353 FlavorHistoryProducer::getClosestJet( Handle<CandidateView> const & pJets,
00354 reco::Candidate const & parton ) const
00355 {
00356 double dr = matchDR_;
00357 CandidateView::const_iterator j = pJets->begin(),
00358 jend = pJets->end();
00359 CandidateView::const_iterator bestJet = pJets->end();
00360 for ( ; j != jend; ++j ) {
00361 double dri = deltaR( parton.p4(), j->p4() );
00362 if ( dri < dr ) {
00363 dr = dri;
00364 bestJet = j;
00365 }
00366 }
00367 return bestJet;
00368 }
00369
00370 #include "FWCore/Framework/interface/MakerMacros.h"
00371
00372 DEFINE_FWK_MODULE( FlavorHistoryProducer );