VHbbDataFormats/interface/JetSubstructureTools.h

#ifndef JETSUBSTRUCTURETOOLS_H
#define JETSUBSTRUCTURETOOLS_H


#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h" 

#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEventAuxInfo.h"
#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEvent.h"

#include <fastjet/ClusterSequence.hh>
#include <fastjet/GhostedAreaSpec.hh>
#include <fastjet/ClusterSequenceArea.hh>
#include "fastjet/tools/Filter.hh"
#include "fastjet/tools/Pruner.hh"
#include "fastjet/tools/MassDropTagger.hh"
#include "fastjet/GhostedAreaSpec.hh"

#include "VHbbAnalysis/VHbbDataFormats/interface/Nsubjettiness.h"

#include <iostream>

class JetSubstructureTools {
    
    // member functions
 public:
    
    // constructor
    JetSubstructureTools( std::vector<float> c_px, std::vector<float> c_py, std::vector<float> c_pz, std::vector<float> c_e, std::vector<float> c_pdgId );    
    ~JetSubstructureTools(){}
    
    fastjet::PseudoJet getPrunedJet(){ return prunedJet_; }
    fastjet::PseudoJet getTrimmedJet(){ return trimmedJet_; }
    fastjet::PseudoJet getFilteredJet(){ return filteredJet_; }    
    float getPrunedJetArea(){ return prunedJetArea_; }
    float getTrimmedJetArea(){ return trimmedJetArea_; }
    float getFilteredJetArea(){ return filteredJetArea_; }
    
    float getTau1(){ return tau1_; }
    float getTau2(){ return tau2_; }
    float getTau3(){ return tau3_; }
    float getTau4(){ return tau4_; }    
    
    fastjet::PseudoJet getPrunedSubJet(int i){ 
        if (i == 0) return prunedSubJet1_;
        else if (i == 1) return prunedSubJet2_;
        else throw cms::Exception("JetSubstructureTools") << "Too many subjets..." << std::endl;
    }

    // data members
 public:
    
    std::vector<fastjet::PseudoJet> FJconstituents_;
    std::vector<float> c_pdgIds_; 
    
    fastjet::PseudoJet prunedJet_;
    fastjet::PseudoJet trimmedJet_;
    fastjet::PseudoJet filteredJet_;
    
    float prunedJetMass_;
    float trimmedJetMass_;
    float filteredJetMass_;
    float ungroomedJetMass_;
    
    float prunedJetArea_;
    float trimmedJetArea_;
    float filteredJetArea_;
    
    float tau1_;
    float tau2_;
    float tau3_;
    float tau4_;
    
    float rcore01_;
    float rcore02_;
    float rcore03_;
    float rcore04_;
    float rcore05_;
    float rcore06_;
    float rcore07_;
    float rcore08_;
    float rcore09_;
    float rcore10_;
    float rcore11_;
    float rcore12_;

    fastjet::PseudoJet prunedSubJet1_;
    fastjet::PseudoJet prunedSubJet2_;    
    
    float QJetVolatility_;
};

    // -------------------------------------------
    // -------------------------------------------
    // -------------------------------------------
    // -------------------------------------------

    // constructor
JetSubstructureTools::JetSubstructureTools( std::vector<float> c_px, std::vector<float> c_py, std::vector<float> c_pz, std::vector<float> c_e, std::vector<float> c_pdgId )   
{
    
    std::cout << "hey y'all!" << std::endl;
    
        // check that they are all the same size
    if ((c_px.size() == c_py.size())&&(c_py.size() == c_pz.size())&&(c_pz.size() == c_e.size())&&(c_e.size() == c_pdgId.size())){
        for (unsigned int i = 0; i < c_px.size(); i++){
            
            FJconstituents_.push_back( fastjet::PseudoJet( c_px[i], c_py[i], c_pz[i], c_e[i] ) );
            c_pdgIds_.push_back( c_pdgId[i] );
            
        }
    }
    else throw cms::Exception("JetSubstructureTools") << "Constituent size mismatch..." << std::endl;
    
        // -------------------------------------------
        // recluster on the fly....
    double mJetRadius = 1.2;
    fastjet::JetDefinition jetDef(fastjet::cambridge_algorithm, mJetRadius);
    
    int activeAreaRepeats = 1;
    double ghostArea = 0.01;
    double ghostEtaMax = 5.0;

    fastjet::GhostedAreaSpec fjActiveArea(ghostEtaMax,activeAreaRepeats,ghostArea);
    fjActiveArea.set_fj2_placement(true);
    fastjet::AreaDefinition fjAreaDefinition( fastjet::active_area_explicit_ghosts, fjActiveArea );

    fastjet::ClusterSequenceArea thisClustering(FJconstituents_, jetDef, fjAreaDefinition);
    std::vector<fastjet::PseudoJet> out_jets = sorted_by_pt(thisClustering.inclusive_jets(50.0));

    fastjet::ClusterSequence thisClustering_basic(FJconstituents_, jetDef);
    std::vector<fastjet::PseudoJet> out_jets_basic = sorted_by_pt(thisClustering_basic.inclusive_jets(50.0));

        // -------------------------------------------
        // define groomers
    fastjet::Filter trimmer( fastjet::Filter(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2), fastjet::SelectorPtFractionMin(0.03)));
    fastjet::Filter filter( fastjet::Filter(fastjet::JetDefinition(fastjet::cambridge_algorithm, 0.3), fastjet::SelectorNHardest(3)));
    fastjet::Pruner pruner(fastjet::cambridge_algorithm, 0.1, 0.5);
    
    std::vector<fastjet::Transformer const *> transformers;
    transformers.push_back(&trimmer);
    transformers.push_back(&filter);
    transformers.push_back(&pruner);

        // define n-subjettiness
    float mNsubjettinessKappa = 1.;
    NsubParameters paraNsub = NsubParameters(mNsubjettinessKappa, mJetRadius);   
    Nsubjettiness routine(nsub_kt_axes, paraNsub);
    
    
    ungroomedJetMass_ = out_jets.at(0).m();
        // -------------------------------------------    
        // compute pruning, trimming, filtering  -------------    
    int transctr = 0;
    for ( std::vector<fastjet::Transformer const *>::const_iterator 
         itransf = transformers.begin(), itransfEnd = transformers.end(); 
         itransf != itransfEnd; ++itransf ) {  
        
        fastjet::PseudoJet transformedJet = out_jets.at(0);
        transformedJet = (**itransf)(transformedJet);
        
        if (transctr == 0){ // trimmed
            trimmedJet_ = transformedJet;
            trimmedJetMass_ = transformedJet.m();            
            trimmedJetArea_ = transformedJet.area();
        }
        else if (transctr == 1){ // filtered
            filteredJet_ = transformedJet;
            filteredJetMass_ = transformedJet.m();            
            filteredJetArea_ = transformedJet.area();
        }
        else if (transctr == 2){ // pruned
            prunedJet_ = transformedJet;
            prunedJetMass_ = transformedJet.m();            
            prunedJetArea_ = transformedJet.area();
            
                //decompose into requested number of subjets:
            if (transformedJet.constituents().size() > 1){
                
                int nsubjetstokeep = 2;
                std::vector<fastjet::PseudoJet> prunedSubjets = transformedJet.associated_cluster_sequence()->exclusive_subjets(transformedJet,nsubjetstokeep);    
                
                prunedSubJet1_ = prunedSubjets.at(0);
                prunedSubJet2_ = prunedSubjets.at(1);
                
            }
        }
        else{ std::cout << "error in number of transformers" << std::endl;}                    
        transctr++;
    }        
    
        // -------------------------------------------    
        // compute n-subjettiness  -------------
    tau1_ = routine.getTau(1, out_jets.at(0).constituents()); 
    tau2_ = routine.getTau(2, out_jets.at(0).constituents());
    tau3_ = routine.getTau(3, out_jets.at(0).constituents());
    tau4_ = routine.getTau(4, out_jets.at(0).constituents());
    
}
    // -------------------------------------------
    // -------------------------------------------
    // -------------------------------------------
    // -------------------------------------------


#endif


Revision:	1.1.2.1
Committed:	Mon Nov 19 20:59:46 2012 UTC (12 years, 5 months ago) by ntran
Content type:	text/plain
Branch:	hbbsubstructDevPostHCP
CVS Tags:	hbbsubstructDev_4, hbbsubstructDev_3, hbbsubstructDev_2
Changes since 1.1:	+214 -0 lines
Log Message:	add files
#	Content
1	#ifndef JETSUBSTRUCTURETOOLS_H
2	#define JETSUBSTRUCTURETOOLS_H
3
4
5	#include "FWCore/ParameterSet/interface/ParameterSet.h"
6	#include "FWCore/Framework/interface/Frameworkfwd.h"
7	#include "FWCore/Framework/interface/EDAnalyzer.h"
8	#include "FWCore/Framework/interface/Event.h"
9
10	#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEventAuxInfo.h"
11	#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEvent.h"
12
13	#include <fastjet/ClusterSequence.hh>
14	#include <fastjet/GhostedAreaSpec.hh>
15	#include <fastjet/ClusterSequenceArea.hh>
16	#include "fastjet/tools/Filter.hh"
17	#include "fastjet/tools/Pruner.hh"
18	#include "fastjet/tools/MassDropTagger.hh"
19	#include "fastjet/GhostedAreaSpec.hh"
20
21	#include "VHbbAnalysis/VHbbDataFormats/interface/Nsubjettiness.h"
22
23	#include <iostream>
24
25	class JetSubstructureTools {
26
27	// member functions
28	public:
29
30	// constructor
31	JetSubstructureTools( std::vector<float> c_px, std::vector<float> c_py, std::vector<float> c_pz, std::vector<float> c_e, std::vector<float> c_pdgId );
32	~JetSubstructureTools(){}
33
34	fastjet::PseudoJet getPrunedJet(){ return prunedJet_; }
35	fastjet::PseudoJet getTrimmedJet(){ return trimmedJet_; }
36	fastjet::PseudoJet getFilteredJet(){ return filteredJet_; }
37	float getPrunedJetArea(){ return prunedJetArea_; }
38	float getTrimmedJetArea(){ return trimmedJetArea_; }
39	float getFilteredJetArea(){ return filteredJetArea_; }
40
41	float getTau1(){ return tau1_; }
42	float getTau2(){ return tau2_; }
43	float getTau3(){ return tau3_; }
44	float getTau4(){ return tau4_; }
45
46	fastjet::PseudoJet getPrunedSubJet(int i){
47	if (i == 0) return prunedSubJet1_;
48	else if (i == 1) return prunedSubJet2_;
49	else throw cms::Exception("JetSubstructureTools") << "Too many subjets..." << std::endl;
50	}
51
52	// data members
53	public:
54
55	std::vector<fastjet::PseudoJet> FJconstituents_;
56	std::vector<float> c_pdgIds_;
57
58	fastjet::PseudoJet prunedJet_;
59	fastjet::PseudoJet trimmedJet_;
60	fastjet::PseudoJet filteredJet_;
61
62	float prunedJetMass_;
63	float trimmedJetMass_;
64	float filteredJetMass_;
65	float ungroomedJetMass_;
66
67	float prunedJetArea_;
68	float trimmedJetArea_;
69	float filteredJetArea_;
70
71	float tau1_;
72	float tau2_;
73	float tau3_;
74	float tau4_;
75
76	float rcore01_;
77	float rcore02_;
78	float rcore03_;
79	float rcore04_;
80	float rcore05_;
81	float rcore06_;
82	float rcore07_;
83	float rcore08_;
84	float rcore09_;
85	float rcore10_;
86	float rcore11_;
87	float rcore12_;
88
89	fastjet::PseudoJet prunedSubJet1_;
90	fastjet::PseudoJet prunedSubJet2_;
91
92	float QJetVolatility_;
93	};
94
95	// -------------------------------------------
96	// -------------------------------------------
97	// -------------------------------------------
98	// -------------------------------------------
99
100	// constructor
101	JetSubstructureTools::JetSubstructureTools( std::vector<float> c_px, std::vector<float> c_py, std::vector<float> c_pz, std::vector<float> c_e, std::vector<float> c_pdgId )
102	{
103
104	std::cout << "hey y'all!" << std::endl;
105
106	// check that they are all the same size
107	if ((c_px.size() == c_py.size())&&(c_py.size() == c_pz.size())&&(c_pz.size() == c_e.size())&&(c_e.size() == c_pdgId.size())){
108	for (unsigned int i = 0; i < c_px.size(); i++){
109
110	FJconstituents_.push_back( fastjet::PseudoJet( c_px[i], c_py[i], c_pz[i], c_e[i] ) );
111	c_pdgIds_.push_back( c_pdgId[i] );
112
113	}
114	}
115	else throw cms::Exception("JetSubstructureTools") << "Constituent size mismatch..." << std::endl;
116
117	// -------------------------------------------
118	// recluster on the fly....
119	double mJetRadius = 1.2;
120	fastjet::JetDefinition jetDef(fastjet::cambridge_algorithm, mJetRadius);
121
122	int activeAreaRepeats = 1;
123	double ghostArea = 0.01;
124	double ghostEtaMax = 5.0;
125
126	fastjet::GhostedAreaSpec fjActiveArea(ghostEtaMax,activeAreaRepeats,ghostArea);
127	fjActiveArea.set_fj2_placement(true);
128	fastjet::AreaDefinition fjAreaDefinition( fastjet::active_area_explicit_ghosts, fjActiveArea );
129
130	fastjet::ClusterSequenceArea thisClustering(FJconstituents_, jetDef, fjAreaDefinition);
131	std::vector<fastjet::PseudoJet> out_jets = sorted_by_pt(thisClustering.inclusive_jets(50.0));
132
133	fastjet::ClusterSequence thisClustering_basic(FJconstituents_, jetDef);
134	std::vector<fastjet::PseudoJet> out_jets_basic = sorted_by_pt(thisClustering_basic.inclusive_jets(50.0));
135
136	// -------------------------------------------
137	// define groomers
138	fastjet::Filter trimmer( fastjet::Filter(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2), fastjet::SelectorPtFractionMin(0.03)));
139	fastjet::Filter filter( fastjet::Filter(fastjet::JetDefinition(fastjet::cambridge_algorithm, 0.3), fastjet::SelectorNHardest(3)));
140	fastjet::Pruner pruner(fastjet::cambridge_algorithm, 0.1, 0.5);
141
142	std::vector<fastjet::Transformer const *> transformers;
143	transformers.push_back(&trimmer);
144	transformers.push_back(&filter);
145	transformers.push_back(&pruner);
146
147	// define n-subjettiness
148	float mNsubjettinessKappa = 1.;
149	NsubParameters paraNsub = NsubParameters(mNsubjettinessKappa, mJetRadius);
150	Nsubjettiness routine(nsub_kt_axes, paraNsub);
151
152
153	ungroomedJetMass_ = out_jets.at(0).m();
154	// -------------------------------------------
155	// compute pruning, trimming, filtering -------------
156	int transctr = 0;
157	for ( std::vector<fastjet::Transformer const *>::const_iterator
158	itransf = transformers.begin(), itransfEnd = transformers.end();
159	itransf != itransfEnd; ++itransf ) {
160
161	fastjet::PseudoJet transformedJet = out_jets.at(0);
162	transformedJet = (**itransf)(transformedJet);
163
164	if (transctr == 0){ // trimmed
165	trimmedJet_ = transformedJet;
166	trimmedJetMass_ = transformedJet.m();
167	trimmedJetArea_ = transformedJet.area();
168	}
169	else if (transctr == 1){ // filtered
170	filteredJet_ = transformedJet;
171	filteredJetMass_ = transformedJet.m();
172	filteredJetArea_ = transformedJet.area();
173	}
174	else if (transctr == 2){ // pruned
175	prunedJet_ = transformedJet;
176	prunedJetMass_ = transformedJet.m();
177	prunedJetArea_ = transformedJet.area();
178
179	//decompose into requested number of subjets:
180	if (transformedJet.constituents().size() > 1){
181
182	int nsubjetstokeep = 2;
183	std::vector<fastjet::PseudoJet> prunedSubjets = transformedJet.associated_cluster_sequence()->exclusive_subjets(transformedJet,nsubjetstokeep);
184
185	prunedSubJet1_ = prunedSubjets.at(0);
186	prunedSubJet2_ = prunedSubjets.at(1);
187
188	}
189	}
190	else{ std::cout << "error in number of transformers" << std::endl;}
191	transctr++;
192	}
193
194	// -------------------------------------------
195	// compute n-subjettiness -------------
196	tau1_ = routine.getTau(1, out_jets.at(0).constituents());
197	tau2_ = routine.getTau(2, out_jets.at(0).constituents());
198	tau3_ = routine.getTau(3, out_jets.at(0).constituents());
199	tau4_ = routine.getTau(4, out_jets.at(0).constituents());
200
201	}
202	// -------------------------------------------
203	// -------------------------------------------
204	// -------------------------------------------
205	// -------------------------------------------
206
207
208	#endif
209
210
211
212
213
214