16 std::vector<float> QjetsVolatility; QjetsVolatility.reserve(jets->size());
23 QjetsVolatility.push_back(-1);
28 vector<fastjet::PseudoJet> allconstits;
30 for (
unsigned k=0;
k < newCand.getJetConstituents().size();
k++){
32 allconstits.push_back( fastjet::PseudoJet( thisParticle->px(), thisParticle->py(), thisParticle->pz(), thisParticle->energy() ) );
35 fastjet::JetDefinition jetDef(fastjet::cambridge_algorithm,
jetRad_);
36 if (
mJetAlgo_==
"AK") jetDef.set_jet_algorithm( fastjet::antikt_algorithm );
37 else if (
mJetAlgo_ ==
"CA") jetDef.set_jet_algorithm( fastjet::cambridge_algorithm );
38 else throw cms::Exception(
"GroomedJetFiller") <<
" unknown jet algorithm " << std::endl;
40 fastjet::ClusterSequence thisClustering_basic(allconstits, jetDef);
41 std::vector<fastjet::PseudoJet> out_jets_basic = sorted_by_pt(thisClustering_basic.inclusive_jets(
cutoff_));
43 if(out_jets_basic.size()!=1){
44 QjetsVolatility.push_back(-1);
51 std::vector<double> qjetmass;
53 vector<fastjet::PseudoJet> constits;
54 unsigned int nqjetconstits = out_jets_basic.at(0).constituents().size();
55 if (nqjetconstits < (
unsigned int)
QjetsPreclustering_) constits = out_jets_basic.at(0).constituents();
56 else constits = out_jets_basic.at(0).associated_cluster_sequence()->exclusive_subjets_up_to(out_jets_basic.at(0),
QjetsPreclustering_);
62 for(
unsigned int ii = 0 ;
ii < (
unsigned int)
ntrial_ ;
ii++){
64 fastjet::ClusterSequence qjet_seq(constits, qjet_def);
65 vector<fastjet::PseudoJet> inclusive_jets2 = sorted_by_pt(qjet_seq.inclusive_jets(
cutoff_));
66 if (inclusive_jets2.size()>0){
67 qjetmass.push_back(inclusive_jets2[0].
m());
72 if (qjetmass.size()<1) {
73 QjetsVolatility.push_back(-1);
77 double mean = std::accumulate( qjetmass.begin( ) , qjetmass.end( ) , 0 ) /qjetmass.size() ;
78 float totalsquared = 0.;
79 for (
unsigned int i = 0;
i < qjetmass.size();
i++){
80 totalsquared += (qjetmass[
i] -
mean)*(qjetmass[
i] - mean) ;
82 float variance =
sqrt( totalsquared/qjetmass.size() );
84 QjetsVolatility.push_back(variance/mean);
89 fillerQJV.insert(jets, QjetsVolatility.begin(), QjetsVolatility.end());
92 iEvent.
put(outQJV,
"QjetsVolatility");
void SetRNEngine(CLHEP::HepRandomEngine *rnEngine)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Base class for all types of Jets.
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
edm::EDGetTokenT< edm::View< reco::Jet > > src_token_
virtual CLHEP::HepRandomEngine & getEngine(StreamID const &) const =0
Use this engine in event methods.
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
StreamID streamID() const