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

  cout << "Main finding routine..." << endl;

  // 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.2
Committed:	Fri May 25 09:31:04 2012 UTC (12 years, 11 months ago) by peiffer
Content type:	text/plain
Branch:	MAIN
Changes since 1.1:	+4 -4 lines
Log Message:	move files to SFrameTools
#	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			cout << "Main finding routine..." << endl;
84
85			// clear the list if it is filled
86			if (JetsOut.size()!=0){
87			JetsOut.clear();
88			}
89
90			unsigned int NumParticles = ParticlesIn.size();
91			if (NumParticles == 0) return;
92
93			// this is the input list for the fastjet finder
94			vector<fastjet::PseudoJet> InputParticles;
95
96			LorentzVector v;
97			// loop over the given particles and fill the input list
98			for (unsigned int ipart=0; ipart<NumParticles; ++ipart){
99			Particle* part = ParticlesIn[ipart];
100			v = part->v4();
101			InputParticles.push_back(fastjet::PseudoJet(v.Px(),v.Py(),v.Pz(),v.E()));
102			}
103
104			// !!!!!!!! run the jet-finder !!!!!!!!!!!!!
105			// delete the old one, if it exists
106			if (fJetFinder){
107			if (!bUseArea) delete fJetFinder;
108			else delete (fastjet::ClusterSequenceArea*) fJetFinder;
109			}
110			if (!bUseArea){
111			fJetFinder = new fastjet::ClusterSequence(InputParticles, *fJetDef);
112			} else {
113			fJetFinder = (fastjet::ClusterSequence) new fastjet::ClusterSequenceArea(InputParticles, fJetDef, *fAreaDef);
114			}
115
116			// get the jets
117			vector<fastjet::PseudoJet> Jets = fJetFinder->inclusive_jets(fPtMin);
118
119			// return if no jets were found
120			if (Jets.size() == 0){
121			return;
122			}
123
124			// sort jets by increasing pt
125			vector<fastjet::PseudoJet> SortedJets = sorted_by_pt(Jets);
126
127			// fill the jets, more sophisticated variables could be calculated (to come)
128			for (unsigned int ijet = 0; ijet < SortedJets.size(); ijet++) {
129			Jet* jet = new Jet();
130			JetsOut.push_back(jet);
131	peiffer	1.2	jet->set_pt ( SortedJets[ijet].pt());
132			jet->set_eta ( SortedJets[ijet].eta());
133			jet->set_phi ( SortedJets[ijet].phi());
134			jet->set_energy ( SortedJets[ijet].E());
135	rkogler	1.1	}
136
137			// get the association of jetparticles to jets
138			fpart_jet_assoc.clear();
139			fpart_jet_assoc = fJetFinder->particle_jet_indices(SortedJets);
140
141			// loop over the input particles and calculate wanted quantities (charged hadron fraction...)
142			// not used at the moment
143			//for (unsigned int ipart=0; ipart<NumParticles; ++ipart){
144			// // particle not inside any jet
145			// if (fpart_jet_assoc[ipart]==-1){
146			// continue;
147			// }
148			// // get the according jet number and store this particle there
149			// int jetnum = fpart_jet_assoc[ipart];
150
151			// if (fpart_jet_assoc[ipart]==IJet){
152			// do something with it
153			// }
154			//}
155
156			return;
157			}
158
159			vector<int> FJet::GetPartJetAssoc()
160			{
161			// returns the jet indices for each particle.
162			// if particle k is in jet j, then the returned vector
163			// has value j at position k
164
165			return fpart_jet_assoc;
166			}
167
168
169			void FJet::Init()
170			{
171			// initialize the jet definition
172			// always performed when a parameter changes
173			if (fJetDef) delete fJetDef;
174			fJetDef = new fastjet::JetDefinition(fJetAlgo, fR, fRecombScheme, fStrategy);
175			}
176
177			void FJet::InitArea()
178			{
179			// initialize the area definition
180			// always performed when an area parameter changed
181
182			if (fAreaDef) delete fAreaDef;
183
184			if (bVoronoiArea){
185			fAreaDef = new fastjet::AreaDefinition(*fVoronoiSpec);
186			} else {
187			fAreaDef = new fastjet::AreaDefinition(fGhostAreaType, *fGhostSpec);
188			}
189
190			}
191
192			void FJet::SetJetAlgorithm(fastjet::JetAlgorithm algo)
193			{
194			// set the jet algorithm,
195			// possibilities: kt, anti-kt, cambridge-aachen
196
197			fJetAlgo = algo;
198			Init();
199			}
200
201
202			void FJet::SetRadius(double radius)
203			{
204			// set the used 'radius'
205			// this is the 1/R term for the kt, anti-kt and Cam/Aachen
206
207			fR = radius;
208			Init();
209
210			}
211
212
213			void FJet::SetRecombScheme(fastjet::RecombinationScheme recom)
214			{
215			// set the recombination scheme,
216			// possibilities: pt, pt2, Et, Et2, E
217
218			fRecombScheme = recom;
219			Init();
220
221			}
222
223			//_________________ getters _______________________
224
225			fastjet::JetAlgorithm FJet::GetJetAlgorithm()
226			{
227			// get the jet algorithm
228			return fJetAlgo;
229			}
230
231
232			fastjet::RecombinationScheme FJet::GetRecombScheme()
233			{
234			// get the recombination scheme, returns unknown if SISCone was used
235			return fRecombScheme;
236			}
237
238			double FJet::GetExclusiveDmerge(int Njet)
239			{
240			// Return the distance parameter dmin from the jet finding algorithm
241			// corresponding to the recombination that went from n+1 to n jets (=Njet)
242			// If the number of particles in the event is <= njets, the function returns 0.
243			// Only call this function after the jet finder has been run
244
245			if (!fJetFinder) return 0;
246			return fJetFinder->exclusive_dmerge(Njet);
247
248			}
249
250
251			//_________________ print information _______________________
252
253			void FJet::Print(const char* usertitle)
254			{
255			// print the settings to the default output
256
257			cout << endl;
258			cout << "+---------------------------------------------------------------" << endl;
259			cout << "\| Settings of the FastJet Interface (" << usertitle << ")" << endl;
260			cout << "+---------------------------------------------------------------" << endl;
261			cout << "\|" << endl << "\| Jet-finder: " << endl;
262			TString des(fJetDef->description());
263			des.ReplaceAll(",", ",\n\| ");
264			cout << "\| " << des << endl;
265			cout << "\| Minimal Pt cut for inclusive jets is set to " << fPtMin << " GeV" << endl;
266
267			cout << "\| Strategy used by FastJet is ";
268			switch(fJetDef->strategy()){
269			case(-4): cout << "N^2 Minimum Heap Tiled (fastest for 500 < N < 10^4)" << endl; break;
270			case(-3): cout << "N^2 Tiled (fastest for 50 < N < 500)" << endl; break;
271			case(-2): cout << "N^2 Poor Tiled (legacy)" << endl; break;
272			case(-1): cout << "N^2 Plain (fastest for N < 50)" << endl; break;
273			case(0): cout << "N^3 Dumb (slowest variant)" << endl; break;
274			case(1): cout << "automatic selection depending on N" << endl; break;
275			case(2): cout << "N ln(N) (fastest for N > 10^4)" << endl; break;
276			case(3): cout << "N ln(N)3pi (legacy)" << endl; break;
277			case(4): cout << "N ln(N)4pi (legacy)" << endl; break;
278			case(12): cout << "N ln(N) Cambridge (exclusively used for cambridge algorithm)" << endl; break;
279			case(13): cout << "N ln(N) Cambridge 2pi2R (exclusively used for cambridge algorithm)" << endl; break;
280			case(14): cout << "N ln(N) Cambridge 4pi (exclusively used for cambridge algorithm)" << endl; break;
281			case(999): cout << "SISCone, which uses its own strategy" << endl; break;
282			default: cout << "unknown strategy" << endl; break;
283			}
284			cout << "\|" << endl;
285			if (!bUseArea){
286			cout << "\| No area calculation was performed." << endl;
287			} else {
288			if (bVoronoiArea){
289			cout << "\| Area calculation was performed with the Voronoi area definition." << endl;
290			cout << "\| "<< fVoronoiSpec->description() << endl;
291			} else {
292			cout << "\| Area calculation was performed with Ghost particles:" << endl;
293			if (fGhostAreaType == fastjet::active_area) cout << "\| Active ghost area was used. " << endl;
294			if (fGhostAreaType == fastjet::active_area_explicit_ghosts) cout << "\| Active ghost area was used, ghost particles were included in output jets." << endl;
295			if (fGhostAreaType == fastjet::one_ghost_passive_area) cout << "\| Passive ghost area was used, clustering the event with one ghost at a time." << endl;
296			if (fGhostAreaType == fastjet::passive_area) cout << "\| Passive ghost area was used, sped up with information specific to the jet-finder algorithm." << endl;
297			TString des(fGhostSpec->description());
298			des.Prepend("\| ");
299			des.ReplaceAll(",", ",\n\| ");
300			cout << des << endl;
301			}
302			}
303			cout << "\|" << endl;
304			cout << "+---------------------------------------------------------------" << endl;
305
306			}