Utilities/src/LargestGenRapidityGap.cc

#include "ForwardAnalysis/Utilities/interface/LargestGenRapidityGap.h"

#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "ForwardAnalysis/Utilities/interface/EtaComparator.h"

#include <vector>
#include <algorithm>

namespace forwardAnalysis {

void LargestGenRapidityGap::operator()(reco::GenParticleCollection const& genParticles, 
                                                    math::XYZTLorentzVector& genGapLow,
                                                    math::XYZTLorentzVector& genGapHigh){
  // Copy and sort gen particles in eta
  std::vector<math::XYZTLorentzVector> genParticlesSort(0);
  reco::GenParticleCollection::const_iterator genpart = genParticles.begin();
  reco::GenParticleCollection::const_iterator genpart_end = genParticles.end();  
  for(; genpart != genpart_end; ++genpart){
        if( genpart->status() != 1 ) continue;

      if((genpart->eta() >= etaEdgeLow_) && (genpart->eta() <= etaEdgeHigh_))
         genParticlesSort.push_back( genpart->p4() );
   }
   std::stable_sort(genParticlesSort.begin(), genParticlesSort.end(), LessByEta<math::XYZTLorentzVector>());

   // Cases: 0, 1 or > 1 particles in selected range
   math::XYZTLorentzVector def_vec(0.,0.,0.,0.);
   if( genParticlesSort.size() == 0 ){
      genGapLow = def_vec; genGapHigh = def_vec;
      return;
   } else if( genParticlesSort.size() == 1 ){
      genGapLow = def_vec;
      genGapHigh = genParticlesSort[0];
      return;
   } else{
      //FIXME; There must be a STL algorithm for this
      double deltaEtaMax = 0.;
      std::vector<math::XYZTLorentzVector>::const_iterator genPartDeltaEtaMax = genParticlesSort.end();
      std::vector<math::XYZTLorentzVector>::const_iterator genpart = genParticlesSort.begin();
      std::vector<math::XYZTLorentzVector>::const_iterator genpart_end = genParticlesSort.end();
      for(; genpart != genpart_end; ++genpart){
         std::vector<math::XYZTLorentzVector>::const_iterator next = genpart + 1;
         double deltaEta = ( next != genpart_end ) ? ( next->eta() - genpart->eta() ) : 0.;
         if( deltaEta > (deltaEtaMax - 0.0001) ){
            deltaEtaMax = deltaEta;
            genPartDeltaEtaMax = genpart;
         } 
      }
      if( genPartDeltaEtaMax != genpart_end ){
         std::vector<math::XYZTLorentzVector>::const_iterator next = genPartDeltaEtaMax + 1;
         if( next != genpart_end ){
            genGapLow = (*genPartDeltaEtaMax);
            genGapHigh = (*next);
         } else{
            genGapLow = def_vec;
            genGapHigh = (*genPartDeltaEtaMax);
         }
      } else{
         genGapLow = def_vec; genGapHigh = def_vec;
         return;
      }
   } 
}

}
Revision:	1.1
Committed:	Wed Aug 24 13:19:04 2011 UTC (13 years, 8 months ago) by antoniov
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V01-01-01, V01-01-00, antoniov-forwardAnalysis-09Jul2012-v1, antoniov-forwardAnalysis-29Jun2012-v1, V01-00-00, antoniov-utilities-11Jun2012-v1, antoniov-forwardAnalysis-Oct072011-v1, sfonseca_10_04_2011, antoniov-forwardAnalysis-Sep182011-v1, antoniov-forwardAnalysis-Sep102011-v1, sfonseca_08_26_2011, forwardAnalysis-Aug232011-v1, HEAD
Error occurred while calculating annotation data.
Log Message:	include LargestGenRapidityGap class
#	Content
1	#include "ForwardAnalysis/Utilities/interface/LargestGenRapidityGap.h"
2
3	#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
4	#include "ForwardAnalysis/Utilities/interface/EtaComparator.h"
5
6	#include <vector>
7	#include <algorithm>
8
9	namespace forwardAnalysis {
10
11	void LargestGenRapidityGap::operator()(reco::GenParticleCollection const& genParticles,
12	math::XYZTLorentzVector& genGapLow,
13	math::XYZTLorentzVector& genGapHigh){
14	// Copy and sort gen particles in eta
15	std::vector<math::XYZTLorentzVector> genParticlesSort(0);
16	reco::GenParticleCollection::const_iterator genpart = genParticles.begin();
17	reco::GenParticleCollection::const_iterator genpart_end = genParticles.end();
18	for(; genpart != genpart_end; ++genpart){
19	if( genpart->status() != 1 ) continue;
20
21	if((genpart->eta() >= etaEdgeLow_) && (genpart->eta() <= etaEdgeHigh_))
22	genParticlesSort.push_back( genpart->p4() );
23	}
24	std::stable_sort(genParticlesSort.begin(), genParticlesSort.end(), LessByEta<math::XYZTLorentzVector>());
25
26	// Cases: 0, 1 or > 1 particles in selected range
27	math::XYZTLorentzVector def_vec(0.,0.,0.,0.);
28	if( genParticlesSort.size() == 0 ){
29	genGapLow = def_vec; genGapHigh = def_vec;
30	return;
31	} else if( genParticlesSort.size() == 1 ){
32	genGapLow = def_vec;
33	genGapHigh = genParticlesSort[0];
34	return;
35	} else{
36	//FIXME; There must be a STL algorithm for this
37	double deltaEtaMax = 0.;
38	std::vector<math::XYZTLorentzVector>::const_iterator genPartDeltaEtaMax = genParticlesSort.end();
39	std::vector<math::XYZTLorentzVector>::const_iterator genpart = genParticlesSort.begin();
40	std::vector<math::XYZTLorentzVector>::const_iterator genpart_end = genParticlesSort.end();
41	for(; genpart != genpart_end; ++genpart){
42	std::vector<math::XYZTLorentzVector>::const_iterator next = genpart + 1;
43	double deltaEta = ( next != genpart_end ) ? ( next->eta() - genpart->eta() ) : 0.;
44	if( deltaEta > (deltaEtaMax - 0.0001) ){
45	deltaEtaMax = deltaEta;
46	genPartDeltaEtaMax = genpart;
47	}
48	}
49	if( genPartDeltaEtaMax != genpart_end ){
50	std::vector<math::XYZTLorentzVector>::const_iterator next = genPartDeltaEtaMax + 1;
51	if( next != genpart_end ){
52	genGapLow = (*genPartDeltaEtaMax);
53	genGapHigh = (*next);
54	} else{
55	genGapLow = def_vec;
56	genGapHigh = (*genPartDeltaEtaMax);
57	}
58	} else{
59	genGapLow = def_vec; genGapHigh = def_vec;
60	return;
61	}
62	}
63	}
64
65	}