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
#	User	Rev	Content
1	rkogler	1.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	peiffer	1.2	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	rkogler	1.1	}
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			}