NtupleWriter/src/JetProps.cxx

#include "UHHAnalysis/NtupleWriter/interface/JetProps.h"

#include <vector>
#include <iostream>

#include "UHHAnalysis/NtupleWriter/interface/Nsubjettiness.h"
#include "UHHAnalysis/NtupleWriter/interface/QjetsPlugin.h"
#include "fastjet/tools/Pruner.hh"

#include "TString.h"
#include "TObject.h"
#include "TMath.h"

using namespace std;

JetProps::JetProps()
{
  // standard constructor
  m_jet = NULL;
  m_pf_candidates = NULL;
  m_JetFinder = NULL;
  m_JetDef = NULL;
}

JetProps::JetProps(TopJet* jet)
{
  // standard constructor
  m_jet = jet;
  m_pf_candidates = NULL;
  m_JetFinder = NULL;
  m_JetDef = NULL;
}

JetProps::~JetProps()
{
  // destructor, clean up
  if (m_JetFinder) delete m_JetFinder;
  if (m_JetDef) delete m_JetDef;
  m_pf_candidates = NULL;
}

double JetProps::GetNsubjettiness(int N, Njettiness::AxesMode mode, double beta, double R0, double Rcutoff)
{
  /// get N-subjettiness
  /// the actual calculation is done by using spartyjet
  /// the second argument is the way the sub-jet axes are found
  /// from NsubjettinessPlugin.hh:
  /// kt_axes              : exclusive kT
  /// ca_axes              : exclusive CA
  /// min_axes             : minimize N-jettiness by iteration (100 passes default)
  /// onepass_{kt,ca}_axes : one-pass minimization seeded by kt or CA (pretty good)

  if (!m_jet){
    return 0.;
  }

  // first: create a fastjet-jet from the TopJet
  std::vector<fastjet::PseudoJet> FJparticles = GetJetConstituents();

  fastjet::PseudoJet fjet = join(FJparticles);
  
  fastjet::Nsubjettiness Nsubj(N, mode, beta, R0, Rcutoff);
  
  double sub=Nsubj.result(fjet);
  
  return sub;

}

double JetProps::GetPrunedNsubjettiness(int N, Njettiness::AxesMode mode, double beta, double R0, double Rcutoff)
{
  /// get N-subjettiness, same as above but for pruned jets
  /// pruning ist first called via GetPrunedJet (see below)

  if (!m_jet){
    return 0.;
  }

  std::vector<fastjet::PseudoJet> jets = GetFastJet(2.0);   // something large to make sure jet is inside radius
  fastjet::PseudoJet pjet = GetPrunedJet(jets[0]);

  // first: create a fastjet-jet from the TopJet
  std::vector<fastjet::PseudoJet> FJparticles = pjet.constituents();
  fastjet::PseudoJet fjet = join(FJparticles);
  
  fastjet::Nsubjettiness Nsubj(N, mode, beta, R0, Rcutoff);
  
  double sub=Nsubj.result(fjet);
  
  return sub;

}

std::vector<fastjet::PseudoJet> JetProps::GetJetConstituents()
{
  // get the jet constituents as Fastjet pseudojets

  std::vector<fastjet::PseudoJet> FJparticles;
  if (m_pf_candidates==NULL){
    cout << "Error, can not get jet constituents without list of PF candidates"<< endl;
    return FJparticles;
  }

  std::vector<unsigned int> jet_parts_ind = m_jet->pfconstituents_indices();

  for(unsigned int ic=0; ic<jet_parts_ind.size(); ++ic){
    PFParticle p = m_pf_candidates->at(jet_parts_ind[ic]);
    fastjet::PseudoJet fj( p.pt()*cos(p.phi()),
                           p.pt()*sin(p.phi()),
                           p.pt()*sinh(p.eta()),
                           p.energy() );
    
    FJparticles.push_back( fj );   
  }
  return FJparticles;
}

vector<fastjet::PseudoJet> JetProps::GetFastJet(double R0, fastjet::JetAlgorithm jet_algo)
{
  /// get fastjet jet from the input particles of the object's jet
  /// could be more than one jet, if the jet definition is different 
  /// than the one being used for the initial jet or the radius is smaller

  std::vector<fastjet::PseudoJet> InputParticles = GetJetConstituents();

  if (m_JetDef){
    delete m_JetDef;
  }
  m_JetDef = new fastjet::JetDefinition(jet_algo, R0);  
  if (m_JetFinder){
    delete m_JetFinder;
  }
  m_JetFinder = new fastjet::ClusterSequence(InputParticles, *m_JetDef);

  vector<fastjet::PseudoJet> Jets = m_JetFinder->inclusive_jets(10.);

  // return if no jets were found
  if (Jets.size() == 0){
    return Jets; 
  }

  // sort jets by increasing pt
  vector<fastjet::PseudoJet> SortedJets = sorted_by_pt(Jets);

  return SortedJets;
  
}

fastjet::PseudoJet JetProps::GetPrunedJet(fastjet::PseudoJet injet)
{
  // calculate the pruned jet

  double zcut = 0.1;   
  double rcut_factor = 0.5;   
  fastjet::Pruner pruner(fastjet::cambridge_algorithm, zcut, rcut_factor);   
  fastjet::PseudoJet pruned_jet = pruner(injet);

  return pruned_jet;
}


double JetProps::GetQjetVolatility(int seed, double R0)
{
    
    fastjet::JetDefinition jetDef(fastjet::cambridge_algorithm, R0);    
    int mQJetsN = 50;
    int mQJetsPreclustering = 30;
    
    std::vector< double > qjetmasses;
    
    std::vector<fastjet::PseudoJet> FJparticles = GetJetConstituents();

    fastjet::ClusterSequence thisClustering_basic(FJparticles, jetDef);
    std::vector<fastjet::PseudoJet> out_jets_basic = sorted_by_pt(thisClustering_basic.inclusive_jets(50.0));
    if (out_jets_basic.size()==0){ 
      return -1.; // no jets found, something wrong
    }

    int j = 0; // the hardest jet
    
    double zcut(0.1), dcut_fctr(0.5), exp_min(0.), exp_max(0.), rigidity(0.1);          
    
    vector<fastjet::PseudoJet> constits;
    unsigned int nqjetconstits = out_jets_basic.at(j).constituents().size();
    if (nqjetconstits < (unsigned int) mQJetsPreclustering) constits = out_jets_basic.at(j).constituents();
    else constits = out_jets_basic.at(j).associated_cluster_sequence()->exclusive_subjets_up_to(out_jets_basic.at(j),mQJetsPreclustering);
    for(unsigned int ii = 0 ; ii < (unsigned int) mQJetsN ; ii++){
        QjetsPlugin qjet_plugin(zcut, dcut_fctr, exp_min, exp_max, rigidity);
        qjet_plugin.SetRandSeed(seed+ii); // new feature in Qjets to set the random seed
        fastjet::JetDefinition qjet_def(&qjet_plugin);
        fastjet::ClusterSequence qjet_seq(constits, qjet_def);
        vector<fastjet::PseudoJet> inclusive_jets2 = sorted_by_pt(qjet_seq.inclusive_jets(50.0));
        
        if (inclusive_jets2.size()>0) { qjetmasses.push_back( inclusive_jets2[0].m() ); }
        
    }
    
        // find RMS of a vector
    double qjetsRMS = FindRMS( qjetmasses );
        // find mean of a vector
    double qjetsMean = FindMean( qjetmasses );
        // compute QJets volatility
    double qjetsVolatility = qjetsRMS/qjetsMean;
    return qjetsVolatility;
}
    // -------------------------------------------
    // -------------------------------------------


double JetProps::FindRMS( std::vector< double > qjetmasses )
{
  
  double total = 0.;
  double ctr = 0.;
  for (unsigned int i = 0; i < qjetmasses.size(); i++){
    total = total + qjetmasses[i];
    ctr++;
  }
  double mean =  total/ctr;
  
  double totalsquared = 0.;
  for (unsigned int i = 0; i < qjetmasses.size(); i++){
    totalsquared += (qjetmasses[i] - mean)*(qjetmasses[i] - mean) ;
  }
  double RMS = sqrt( totalsquared/ctr );
  return RMS;
}

double JetProps::FindMean( std::vector< double > qjetmasses )
{
  double total = 0.;
  double ctr = 0.;
  for (unsigned int i = 0; i < qjetmasses.size(); i++){
    total = total + qjetmasses[i];
    ctr++;
  }
  return total/ctr;
}
Revision:	1.2
Committed:	Wed Jun 19 22:09:57 2013 UTC (11 years, 10 months ago) by rkogler
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.1:	+0 -9 lines
Log Message:	removed output statement
#	User	Rev	Content
1	rkogler	1.1	#include "UHHAnalysis/NtupleWriter/interface/JetProps.h"
2
3			#include <vector>
4			#include <iostream>
5
6			#include "UHHAnalysis/NtupleWriter/interface/Nsubjettiness.h"
7			#include "UHHAnalysis/NtupleWriter/interface/QjetsPlugin.h"
8			#include "fastjet/tools/Pruner.hh"
9
10			#include "TString.h"
11			#include "TObject.h"
12			#include "TMath.h"
13
14			using namespace std;
15
16			JetProps::JetProps()
17			{
18			// standard constructor
19			m_jet = NULL;
20			m_pf_candidates = NULL;
21			m_JetFinder = NULL;
22			m_JetDef = NULL;
23			}
24
25			JetProps::JetProps(TopJet* jet)
26			{
27			// standard constructor
28			m_jet = jet;
29			m_pf_candidates = NULL;
30			m_JetFinder = NULL;
31			m_JetDef = NULL;
32			}
33
34			JetProps::~JetProps()
35			{
36			// destructor, clean up
37			if (m_JetFinder) delete m_JetFinder;
38			if (m_JetDef) delete m_JetDef;
39			m_pf_candidates = NULL;
40			}
41
42			double JetProps::GetNsubjettiness(int N, Njettiness::AxesMode mode, double beta, double R0, double Rcutoff)
43			{
44			/// get N-subjettiness
45			/// the actual calculation is done by using spartyjet
46			/// the second argument is the way the sub-jet axes are found
47			/// from NsubjettinessPlugin.hh:
48			/// kt_axes : exclusive kT
49			/// ca_axes : exclusive CA
50			/// min_axes : minimize N-jettiness by iteration (100 passes default)
51			/// onepass_{kt,ca}_axes : one-pass minimization seeded by kt or CA (pretty good)
52
53			if (!m_jet){
54			return 0.;
55			}
56
57			// first: create a fastjet-jet from the TopJet
58			std::vector<fastjet::PseudoJet> FJparticles = GetJetConstituents();
59
60			fastjet::PseudoJet fjet = join(FJparticles);
61
62			fastjet::Nsubjettiness Nsubj(N, mode, beta, R0, Rcutoff);
63
64			double sub=Nsubj.result(fjet);
65
66			return sub;
67
68			}
69
70			double JetProps::GetPrunedNsubjettiness(int N, Njettiness::AxesMode mode, double beta, double R0, double Rcutoff)
71			{
72			/// get N-subjettiness, same as above but for pruned jets
73			/// pruning ist first called via GetPrunedJet (see below)
74
75			if (!m_jet){
76			return 0.;
77			}
78
79			std::vector<fastjet::PseudoJet> jets = GetFastJet(2.0); // something large to make sure jet is inside radius
80			fastjet::PseudoJet pjet = GetPrunedJet(jets[0]);
81
82			// first: create a fastjet-jet from the TopJet
83			std::vector<fastjet::PseudoJet> FJparticles = pjet.constituents();
84			fastjet::PseudoJet fjet = join(FJparticles);
85
86			fastjet::Nsubjettiness Nsubj(N, mode, beta, R0, Rcutoff);
87
88			double sub=Nsubj.result(fjet);
89
90			return sub;
91
92			}
93
94			std::vector<fastjet::PseudoJet> JetProps::GetJetConstituents()
95			{
96			// get the jet constituents as Fastjet pseudojets
97
98			std::vector<fastjet::PseudoJet> FJparticles;
99			if (m_pf_candidates==NULL){
100			cout << "Error, can not get jet constituents without list of PF candidates"<< endl;
101			return FJparticles;
102			}
103
104			std::vector<unsigned int> jet_parts_ind = m_jet->pfconstituents_indices();
105
106			for(unsigned int ic=0; ic<jet_parts_ind.size(); ++ic){
107			PFParticle p = m_pf_candidates->at(jet_parts_ind[ic]);
108			fastjet::PseudoJet fj( p.pt()*cos(p.phi()),
109			p.pt()*sin(p.phi()),
110			p.pt()*sinh(p.eta()),
111			p.energy() );
112
113			FJparticles.push_back( fj );
114			}
115			return FJparticles;
116			}
117
118			vector<fastjet::PseudoJet> JetProps::GetFastJet(double R0, fastjet::JetAlgorithm jet_algo)
119			{
120			/// get fastjet jet from the input particles of the object's jet
121			/// could be more than one jet, if the jet definition is different
122			/// than the one being used for the initial jet or the radius is smaller
123
124			std::vector<fastjet::PseudoJet> InputParticles = GetJetConstituents();
125
126			if (m_JetDef){
127			delete m_JetDef;
128			}
129			m_JetDef = new fastjet::JetDefinition(jet_algo, R0);
130			if (m_JetFinder){
131			delete m_JetFinder;
132			}
133			m_JetFinder = new fastjet::ClusterSequence(InputParticles, *m_JetDef);
134
135			vector<fastjet::PseudoJet> Jets = m_JetFinder->inclusive_jets(10.);
136
137			// return if no jets were found
138			if (Jets.size() == 0){
139			return Jets;
140			}
141
142			// sort jets by increasing pt
143			vector<fastjet::PseudoJet> SortedJets = sorted_by_pt(Jets);
144
145			return SortedJets;
146
147			}
148
149			fastjet::PseudoJet JetProps::GetPrunedJet(fastjet::PseudoJet injet)
150			{
151			// calculate the pruned jet
152
153			double zcut = 0.1;
154			double rcut_factor = 0.5;
155			fastjet::Pruner pruner(fastjet::cambridge_algorithm, zcut, rcut_factor);
156			fastjet::PseudoJet pruned_jet = pruner(injet);
157
158			return pruned_jet;
159			}
160
161
162
163			double JetProps::GetQjetVolatility(int seed, double R0)
164			{
165
166			fastjet::JetDefinition jetDef(fastjet::cambridge_algorithm, R0);
167			int mQJetsN = 50;
168			int mQJetsPreclustering = 30;
169
170			std::vector< double > qjetmasses;
171
172			std::vector<fastjet::PseudoJet> FJparticles = GetJetConstituents();
173
174			fastjet::ClusterSequence thisClustering_basic(FJparticles, jetDef);
175			std::vector<fastjet::PseudoJet> out_jets_basic = sorted_by_pt(thisClustering_basic.inclusive_jets(50.0));
176			if (out_jets_basic.size()==0){
177			return -1.; // no jets found, something wrong
178			}
179
180			int j = 0; // the hardest jet
181
182			double zcut(0.1), dcut_fctr(0.5), exp_min(0.), exp_max(0.), rigidity(0.1);
183
184			vector<fastjet::PseudoJet> constits;
185			unsigned int nqjetconstits = out_jets_basic.at(j).constituents().size();
186			if (nqjetconstits < (unsigned int) mQJetsPreclustering) constits = out_jets_basic.at(j).constituents();
187			else constits = out_jets_basic.at(j).associated_cluster_sequence()->exclusive_subjets_up_to(out_jets_basic.at(j),mQJetsPreclustering);
188			for(unsigned int ii = 0 ; ii < (unsigned int) mQJetsN ; ii++){
189			QjetsPlugin qjet_plugin(zcut, dcut_fctr, exp_min, exp_max, rigidity);
190			qjet_plugin.SetRandSeed(seed+ii); // new feature in Qjets to set the random seed
191			fastjet::JetDefinition qjet_def(&qjet_plugin);
192			fastjet::ClusterSequence qjet_seq(constits, qjet_def);
193			vector<fastjet::PseudoJet> inclusive_jets2 = sorted_by_pt(qjet_seq.inclusive_jets(50.0));
194
195			if (inclusive_jets2.size()>0) { qjetmasses.push_back( inclusive_jets2[0].m() ); }
196
197			}
198
199			// find RMS of a vector
200			double qjetsRMS = FindRMS( qjetmasses );
201			// find mean of a vector
202			double qjetsMean = FindMean( qjetmasses );
203			// compute QJets volatility
204			double qjetsVolatility = qjetsRMS/qjetsMean;
205			return qjetsVolatility;
206			}
207			// -------------------------------------------
208			// -------------------------------------------
209
210
211
212			double JetProps::FindRMS( std::vector< double > qjetmasses )
213			{
214
215			double total = 0.;
216			double ctr = 0.;
217			for (unsigned int i = 0; i < qjetmasses.size(); i++){
218			total = total + qjetmasses[i];
219			ctr++;
220			}
221			double mean = total/ctr;
222
223			double totalsquared = 0.;
224			for (unsigned int i = 0; i < qjetmasses.size(); i++){
225			totalsquared += (qjetmasses[i] - mean)*(qjetmasses[i] - mean) ;
226			}
227			double RMS = sqrt( totalsquared/ctr );
228			return RMS;
229			}
230
231			double JetProps::FindMean( std::vector< double > qjetmasses )
232			{
233			double total = 0.;
234			double ctr = 0.;
235			for (unsigned int i = 0; i < qjetmasses.size(); i++){
236			total = total + qjetmasses[i];
237			ctr++;
238			}
239			return total/ctr;
240			}