SFrameTools/src/FJet.cxx

#include "../include/FJet.h"
#include "../include/Objects.h"

#include <vector>
#include "Math/LorentzVector.h"
#include "Math/PtEtaPhiE4D.h"
#include "TString.h"
#include "TObject.h"
#include "TMath.h"

typedef ROOT::Math::LorentzVector< ROOT::Math::PtEtaPhiE4D< Double32_t > > LorentzVector;
using namespace std;

FJet::FJet()
{
  // standard constructor, use the following as defaults:
  // jet algorithm = Cam/Aachen
  // radius = 0.8
  // Pt_min = 3.0 GeV
  // recombination scheme = E-scheme (four-vector addition)
  // clustering strategy = best 

  fJetAlgo = fastjet::cambridge_algorithm;
  //fJetAlgo = fastjet::kt_algorithm;
  fR = 0.8;
  fPtMin = 10.0;
  fRecombScheme = fastjet::E_scheme; 
  fStrategy = fastjet::Best;
  fJetDef = new fastjet::JetDefinition(fJetAlgo, fR, fRecombScheme, fStrategy);
   
  // the jet finder, initialise it to NULL
  fJetFinder = NULL;
  
  // default settings for the area treatment
  bUseArea = false;
  fAreaDef = NULL;
  fGhostSpec = new fastjet::GhostedAreaSpec(4.1, 1, 0.02, 1.0, 0.1, 1e-100);
  fVoronoiSpec = new fastjet::VoronoiAreaSpec();
  fGhostAreaType = fastjet::active_area;
  bVoronoiArea = false;

}

FJet::~FJet()
{
  // destructor, clean up
  if (fJetDef) delete fJetDef;
  if (fAreaDef) delete fAreaDef;
  if (fGhostSpec) delete fGhostSpec;
  if (fVoronoiSpec) delete fVoronoiSpec;
  fJetDef=NULL;
  fAreaDef=NULL;
  fGhostSpec=NULL;
  fVoronoiSpec=NULL;

}

void FJet::PrepareInput(vector<GenParticle> genparts, vector<Particle*>& out)
{
  // fill the input array for the jetfinder with 
  // pointers to objects of type Particle, a memcopy is used

  if (out.size()!=0){
    out.clear();
  }

  for (unsigned int i=0; i<genparts.size(); ++i){
    GenParticle* gen = new GenParticle();
    memcpy(gen, &genparts[i], sizeof(genparts[i]));
    Particle* part = (Particle*)gen;
    out.push_back(part);
  }
  return;

}

void FJet::FindJets(const vector<Particle*> ParticlesIn, vector<Jet*>& JetsOut)
{

  // the main routine:
  // jet finding using FastJet

  // clear the list if it is filled
  if (JetsOut.size()!=0){
    JetsOut.clear();
  }
    
  unsigned int NumParticles = ParticlesIn.size();
  if (NumParticles == 0) return;  

  // this is the input list for the fastjet finder
  vector<fastjet::PseudoJet> InputParticles;

  LorentzVector v;
  // loop over the given particles and fill the input list
  for (unsigned int ipart=0; ipart<NumParticles; ++ipart){
    Particle* part = ParticlesIn[ipart];
    v = part->v4();
    InputParticles.push_back(fastjet::PseudoJet(v.Px(),v.Py(),v.Pz(),v.E()));     
  }

  // !!!!!!!!  run the jet-finder !!!!!!!!!!!!!
  // delete the old one, if it exists
  if (fJetFinder){
    if (!bUseArea) delete fJetFinder;
    else delete (fastjet::ClusterSequenceArea*) fJetFinder; 
  }
  if (!bUseArea){
    fJetFinder = new fastjet::ClusterSequence(InputParticles, *fJetDef);
  } else {
    fJetFinder = (fastjet::ClusterSequence*) new fastjet::ClusterSequenceArea(InputParticles, *fJetDef, *fAreaDef);
  }

  // get the jets
  vector<fastjet::PseudoJet> Jets = fJetFinder->inclusive_jets(fPtMin);

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

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

  // fill the jets, more sophisticated variables could be calculated (to come)
  for (unsigned int ijet = 0; ijet < SortedJets.size(); ijet++) { 
    Jet* jet = new Jet();
    JetsOut.push_back(jet);
    jet->set_pt ( SortedJets[ijet].pt());
    jet->set_eta ( SortedJets[ijet].eta());
    jet->set_phi ( SortedJets[ijet].phi());
    jet->set_energy ( SortedJets[ijet].E());
  }

  // get the association of jetparticles to jets
  fpart_jet_assoc.clear();
  fpart_jet_assoc = fJetFinder->particle_jet_indices(SortedJets);

  // loop over the input particles and calculate wanted quantities (charged hadron fraction...)
  // not used at the moment
  //for (unsigned int ipart=0; ipart<NumParticles; ++ipart){
  //     // particle not inside any jet
  //     if (fpart_jet_assoc[ipart]==-1){
  //       continue;
  //     }    
  //     // get the according jet number and store this particle there
  //     int jetnum = fpart_jet_assoc[ipart];

  //  if (fpart_jet_assoc[ipart]==IJet){
    // do something with it
  //  }
  //}
  
  return;
}

vector<int> FJet::GetPartJetAssoc()
{
  // returns the jet indices for each particle. 
  // if particle k is in jet j, then the returned vector 
  // has value j at position k

  return fpart_jet_assoc;  
}


void FJet::Init()
{
  // initialize the jet definition 
  // always performed when a parameter changes
  if (fJetDef) delete fJetDef;
  fJetDef = new fastjet::JetDefinition(fJetAlgo, fR, fRecombScheme, fStrategy);
}

void FJet::InitArea()
{
  // initialize the area definition
  // always performed when an area parameter changed
  
  if (fAreaDef) delete fAreaDef;

  if (bVoronoiArea){
    fAreaDef = new fastjet::AreaDefinition(*fVoronoiSpec);
  } else {
    fAreaDef = new fastjet::AreaDefinition(fGhostAreaType, *fGhostSpec);
  }

}

void FJet::SetJetAlgorithm(fastjet::JetAlgorithm algo)
{
  // set the jet algorithm, 
  // possibilities: kt, anti-kt, cambridge-aachen
  
  fJetAlgo = algo;
  Init();   
}


void FJet::SetRadius(double radius)
{
  // set the used 'radius' 
  // this is the 1/R term for the kt, anti-kt and Cam/Aachen

  fR = radius; 
  Init(); 

}


void FJet::SetRecombScheme(fastjet::RecombinationScheme recom)
{
  // set the recombination scheme, 
  // possibilities: pt, pt2, Et, Et2, E

  fRecombScheme = recom; 
  Init(); 

}

//_________________  getters _______________________

fastjet::JetAlgorithm FJet::GetJetAlgorithm()
{
  // get the jet algorithm
  return fJetAlgo;
}


fastjet::RecombinationScheme FJet::GetRecombScheme()
{
  // get the recombination scheme, returns unknown if SISCone was used
  return fRecombScheme;
}

double FJet::GetExclusiveDmerge(int Njet)
{
  // Return the distance parameter dmin from the jet finding algorithm 
  // corresponding to the recombination that went from n+1 to n jets (=Njet)
  // If the number of particles in the event is <= njets, the function returns 0. 
  // Only call this function after the jet finder has been run

  if (!fJetFinder) return 0;
  return fJetFinder->exclusive_dmerge(Njet);

}


//_________________  print information _______________________

void FJet::Print(const char* usertitle)
{
  // print the settings to the default output

  cout << endl;
  cout << "+---------------------------------------------------------------" << endl;
  cout << "| Settings of the FastJet Interface (" << usertitle << ")" << endl;
  cout << "+---------------------------------------------------------------" << endl; 
  cout << "|" << endl << "| Jet-finder: " << endl;
  TString des(fJetDef->description());
  des.ReplaceAll(",", ",\n| ");
  cout << "|  " << des << endl;
  cout << "|  Minimal Pt cut for inclusive jets is set to " << fPtMin << " GeV" << endl;

  cout << "|  Strategy used by FastJet is ";
  switch(fJetDef->strategy()){
  case(-4): cout << "N^2 Minimum Heap Tiled (fastest for 500 < N < 10^4)" << endl; break;
  case(-3): cout << "N^2 Tiled (fastest for 50 < N < 500)" << endl; break;
  case(-2): cout << "N^2 Poor Tiled (legacy)" << endl; break;
  case(-1): cout << "N^2 Plain (fastest for N < 50)" << endl; break;
  case(0): cout << "N^3 Dumb (slowest variant)" << endl; break;
  case(1): cout << "automatic selection depending on N" << endl; break;
  case(2): cout << "N ln(N) (fastest for N > 10^4)" << endl; break;
  case(3): cout << "N ln(N)3pi (legacy)" << endl; break;
  case(4): cout << "N ln(N)4pi (legacy)" << endl; break;
  case(12): cout << "N ln(N) Cambridge (exclusively used for cambridge algorithm)" << endl; break;
  case(13): cout << "N ln(N) Cambridge 2pi2R (exclusively used for cambridge algorithm)" << endl; break;
  case(14): cout << "N ln(N) Cambridge 4pi (exclusively used for cambridge algorithm)" << endl; break;
  case(999): cout << "SISCone, which uses its own strategy" << endl; break;
  default: cout << "unknown strategy" << endl; break;
  }
  cout << "|" << endl;
  if (!bUseArea){
    cout << "| No area calculation was performed." << endl; 
  } else {
    if (bVoronoiArea){
      cout << "| Area calculation was performed with the Voronoi area definition." << endl; 
      cout << "|  "<< fVoronoiSpec->description() << endl;
    } else {
      cout << "| Area calculation was performed with Ghost particles:" << endl;
      if (fGhostAreaType == fastjet::active_area) cout << "|  Active ghost area was used. " << endl;
      if (fGhostAreaType == fastjet::active_area_explicit_ghosts) cout << "|  Active ghost area was used, ghost particles were included in output jets." << endl;
      if (fGhostAreaType == fastjet::one_ghost_passive_area) cout << "|  Passive ghost area was used, clustering the event with one ghost at a time." << endl;
      if (fGhostAreaType == fastjet::passive_area) cout << "|  Passive ghost area was used, sped up with information specific to the jet-finder algorithm." << endl;
      TString des(fGhostSpec->description());
      des.Prepend("|  ");
      des.ReplaceAll(",", ",\n| ");
      cout << des << endl;
    }
  }
  cout << "|" << endl;
  cout << "+---------------------------------------------------------------" << endl; 

}
Revision:	1.3
Committed:	Wed Jun 6 08:49:57 2012 UTC (12 years, 11 months ago) by rkogler
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Makefile, v1-00, Feb-15-2013-v1, Feb-14-2013, Feb-07-2013-v1, Jan-17-2013-v2, Jan-17-2013-v1, Jan-16-2012-v1, Jan-09-2012-v2, Jan-09-2012-v1, Dec-26-2012-v1, Dec-20-2012-v1, Dec-17-2012-v1, Nov-30-2012-v2, Nov-30-2012-v1, HEAD
Changes since 1.2:	+0 -2 lines
Log Message:	printout removed
#	Content
1	#include "../include/FJet.h"
2	#include "../include/Objects.h"
3
4	#include <vector>
5	#include "Math/LorentzVector.h"
6	#include "Math/PtEtaPhiE4D.h"
7	#include "TString.h"
8	#include "TObject.h"
9	#include "TMath.h"
10
11	typedef ROOT::Math::LorentzVector< ROOT::Math::PtEtaPhiE4D< Double32_t > > LorentzVector;
12	using namespace std;
13
14	FJet::FJet()
15	{
16	// standard constructor, use the following as defaults:
17	// jet algorithm = Cam/Aachen
18	// radius = 0.8
19	// Pt_min = 3.0 GeV
20	// recombination scheme = E-scheme (four-vector addition)
21	// clustering strategy = best
22
23	fJetAlgo = fastjet::cambridge_algorithm;
24	//fJetAlgo = fastjet::kt_algorithm;
25	fR = 0.8;
26	fPtMin = 10.0;
27	fRecombScheme = fastjet::E_scheme;
28	fStrategy = fastjet::Best;
29	fJetDef = new fastjet::JetDefinition(fJetAlgo, fR, fRecombScheme, fStrategy);
30
31	// the jet finder, initialise it to NULL
32	fJetFinder = NULL;
33
34	// default settings for the area treatment
35	bUseArea = false;
36	fAreaDef = NULL;
37	fGhostSpec = new fastjet::GhostedAreaSpec(4.1, 1, 0.02, 1.0, 0.1, 1e-100);
38	fVoronoiSpec = new fastjet::VoronoiAreaSpec();
39	fGhostAreaType = fastjet::active_area;
40	bVoronoiArea = false;
41
42	}
43
44	FJet::~FJet()
45	{
46	// destructor, clean up
47	if (fJetDef) delete fJetDef;
48	if (fAreaDef) delete fAreaDef;
49	if (fGhostSpec) delete fGhostSpec;
50	if (fVoronoiSpec) delete fVoronoiSpec;
51	fJetDef=NULL;
52	fAreaDef=NULL;
53	fGhostSpec=NULL;
54	fVoronoiSpec=NULL;
55
56	}
57
58	void FJet::PrepareInput(vector<GenParticle> genparts, vector<Particle*>& out)
59	{
60	// fill the input array for the jetfinder with
61	// pointers to objects of type Particle, a memcopy is used
62
63	if (out.size()!=0){
64	out.clear();
65	}
66
67	for (unsigned int i=0; i<genparts.size(); ++i){
68	GenParticle* gen = new GenParticle();
69	memcpy(gen, &genparts[i], sizeof(genparts[i]));
70	Particle* part = (Particle*)gen;
71	out.push_back(part);
72	}
73	return;
74
75	}
76
77	void FJet::FindJets(const vector<Particle> ParticlesIn, vector<Jet>& JetsOut)
78	{
79
80	// the main routine:
81	// jet finding using FastJet
82
83	// clear the list if it is filled
84	if (JetsOut.size()!=0){
85	JetsOut.clear();
86	}
87
88	unsigned int NumParticles = ParticlesIn.size();
89	if (NumParticles == 0) return;
90
91	// this is the input list for the fastjet finder
92	vector<fastjet::PseudoJet> InputParticles;
93
94	LorentzVector v;
95	// loop over the given particles and fill the input list
96	for (unsigned int ipart=0; ipart<NumParticles; ++ipart){
97	Particle* part = ParticlesIn[ipart];
98	v = part->v4();
99	InputParticles.push_back(fastjet::PseudoJet(v.Px(),v.Py(),v.Pz(),v.E()));
100	}
101
102	// !!!!!!!! run the jet-finder !!!!!!!!!!!!!
103	// delete the old one, if it exists
104	if (fJetFinder){
105	if (!bUseArea) delete fJetFinder;
106	else delete (fastjet::ClusterSequenceArea*) fJetFinder;
107	}
108	if (!bUseArea){
109	fJetFinder = new fastjet::ClusterSequence(InputParticles, *fJetDef);
110	} else {
111	fJetFinder = (fastjet::ClusterSequence) new fastjet::ClusterSequenceArea(InputParticles, fJetDef, *fAreaDef);
112	}
113
114	// get the jets
115	vector<fastjet::PseudoJet> Jets = fJetFinder->inclusive_jets(fPtMin);
116
117	// return if no jets were found
118	if (Jets.size() == 0){
119	return;
120	}
121
122	// sort jets by increasing pt
123	vector<fastjet::PseudoJet> SortedJets = sorted_by_pt(Jets);
124
125	// fill the jets, more sophisticated variables could be calculated (to come)
126	for (unsigned int ijet = 0; ijet < SortedJets.size(); ijet++) {
127	Jet* jet = new Jet();
128	JetsOut.push_back(jet);
129	jet->set_pt ( SortedJets[ijet].pt());
130	jet->set_eta ( SortedJets[ijet].eta());
131	jet->set_phi ( SortedJets[ijet].phi());
132	jet->set_energy ( SortedJets[ijet].E());
133	}
134
135	// get the association of jetparticles to jets
136	fpart_jet_assoc.clear();
137	fpart_jet_assoc = fJetFinder->particle_jet_indices(SortedJets);
138
139	// loop over the input particles and calculate wanted quantities (charged hadron fraction...)
140	// not used at the moment
141	//for (unsigned int ipart=0; ipart<NumParticles; ++ipart){
142	// // particle not inside any jet
143	// if (fpart_jet_assoc[ipart]==-1){
144	// continue;
145	// }
146	// // get the according jet number and store this particle there
147	// int jetnum = fpart_jet_assoc[ipart];
148
149	// if (fpart_jet_assoc[ipart]==IJet){
150	// do something with it
151	// }
152	//}
153
154	return;
155	}
156
157	vector<int> FJet::GetPartJetAssoc()
158	{
159	// returns the jet indices for each particle.
160	// if particle k is in jet j, then the returned vector
161	// has value j at position k
162
163	return fpart_jet_assoc;
164	}
165
166
167	void FJet::Init()
168	{
169	// initialize the jet definition
170	// always performed when a parameter changes
171	if (fJetDef) delete fJetDef;
172	fJetDef = new fastjet::JetDefinition(fJetAlgo, fR, fRecombScheme, fStrategy);
173	}
174
175	void FJet::InitArea()
176	{
177	// initialize the area definition
178	// always performed when an area parameter changed
179
180	if (fAreaDef) delete fAreaDef;
181
182	if (bVoronoiArea){
183	fAreaDef = new fastjet::AreaDefinition(*fVoronoiSpec);
184	} else {
185	fAreaDef = new fastjet::AreaDefinition(fGhostAreaType, *fGhostSpec);
186	}
187
188	}
189
190	void FJet::SetJetAlgorithm(fastjet::JetAlgorithm algo)
191	{
192	// set the jet algorithm,
193	// possibilities: kt, anti-kt, cambridge-aachen
194
195	fJetAlgo = algo;
196	Init();
197	}
198
199
200	void FJet::SetRadius(double radius)
201	{
202	// set the used 'radius'
203	// this is the 1/R term for the kt, anti-kt and Cam/Aachen
204
205	fR = radius;
206	Init();
207
208	}
209
210
211	void FJet::SetRecombScheme(fastjet::RecombinationScheme recom)
212	{
213	// set the recombination scheme,
214	// possibilities: pt, pt2, Et, Et2, E
215
216	fRecombScheme = recom;
217	Init();
218
219	}
220
221	//_________________ getters _______________________
222
223	fastjet::JetAlgorithm FJet::GetJetAlgorithm()
224	{
225	// get the jet algorithm
226	return fJetAlgo;
227	}
228
229
230	fastjet::RecombinationScheme FJet::GetRecombScheme()
231	{
232	// get the recombination scheme, returns unknown if SISCone was used
233	return fRecombScheme;
234	}
235
236	double FJet::GetExclusiveDmerge(int Njet)
237	{
238	// Return the distance parameter dmin from the jet finding algorithm
239	// corresponding to the recombination that went from n+1 to n jets (=Njet)
240	// If the number of particles in the event is <= njets, the function returns 0.
241	// Only call this function after the jet finder has been run
242
243	if (!fJetFinder) return 0;
244	return fJetFinder->exclusive_dmerge(Njet);
245
246	}
247
248
249	//_________________ print information _______________________
250
251	void FJet::Print(const char* usertitle)
252	{
253	// print the settings to the default output
254
255	cout << endl;
256	cout << "+---------------------------------------------------------------" << endl;
257	cout << "\| Settings of the FastJet Interface (" << usertitle << ")" << endl;
258	cout << "+---------------------------------------------------------------" << endl;
259	cout << "\|" << endl << "\| Jet-finder: " << endl;
260	TString des(fJetDef->description());
261	des.ReplaceAll(",", ",\n\| ");
262	cout << "\| " << des << endl;
263	cout << "\| Minimal Pt cut for inclusive jets is set to " << fPtMin << " GeV" << endl;
264
265	cout << "\| Strategy used by FastJet is ";
266	switch(fJetDef->strategy()){
267	case(-4): cout << "N^2 Minimum Heap Tiled (fastest for 500 < N < 10^4)" << endl; break;
268	case(-3): cout << "N^2 Tiled (fastest for 50 < N < 500)" << endl; break;
269	case(-2): cout << "N^2 Poor Tiled (legacy)" << endl; break;
270	case(-1): cout << "N^2 Plain (fastest for N < 50)" << endl; break;
271	case(0): cout << "N^3 Dumb (slowest variant)" << endl; break;
272	case(1): cout << "automatic selection depending on N" << endl; break;
273	case(2): cout << "N ln(N) (fastest for N > 10^4)" << endl; break;
274	case(3): cout << "N ln(N)3pi (legacy)" << endl; break;
275	case(4): cout << "N ln(N)4pi (legacy)" << endl; break;
276	case(12): cout << "N ln(N) Cambridge (exclusively used for cambridge algorithm)" << endl; break;
277	case(13): cout << "N ln(N) Cambridge 2pi2R (exclusively used for cambridge algorithm)" << endl; break;
278	case(14): cout << "N ln(N) Cambridge 4pi (exclusively used for cambridge algorithm)" << endl; break;
279	case(999): cout << "SISCone, which uses its own strategy" << endl; break;
280	default: cout << "unknown strategy" << endl; break;
281	}
282	cout << "\|" << endl;
283	if (!bUseArea){
284	cout << "\| No area calculation was performed." << endl;
285	} else {
286	if (bVoronoiArea){
287	cout << "\| Area calculation was performed with the Voronoi area definition." << endl;
288	cout << "\| "<< fVoronoiSpec->description() << endl;
289	} else {
290	cout << "\| Area calculation was performed with Ghost particles:" << endl;
291	if (fGhostAreaType == fastjet::active_area) cout << "\| Active ghost area was used. " << endl;
292	if (fGhostAreaType == fastjet::active_area_explicit_ghosts) cout << "\| Active ghost area was used, ghost particles were included in output jets." << endl;
293	if (fGhostAreaType == fastjet::one_ghost_passive_area) cout << "\| Passive ghost area was used, clustering the event with one ghost at a time." << endl;
294	if (fGhostAreaType == fastjet::passive_area) cout << "\| Passive ghost area was used, sped up with information specific to the jet-finder algorithm." << endl;
295	TString des(fGhostSpec->description());
296	des.Prepend("\| ");
297	des.ReplaceAll(",", ",\n\| ");
298	cout << des << endl;
299	}
300	}
301	cout << "\|" << endl;
302	cout << "+---------------------------------------------------------------" << endl;
303
304	}