ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
hijfst.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file hijfst.cc
1 #include <fastjet/ClusterSequence.hh>
2 
3 #include <algorithm>
4 #include <map>
5 #include <string>
6 #include <vector>
7 
8 using namespace std;
9 using namespace fastjet;
10 
11 bool enablep = false;
12 struct algo_info
13 {
14  JetAlgorithm algorithm;
15  double R;
16  int PID;
17  double EtaMin;
18  double EtaMax;
19  double EtMin;
20 };
21 vector<algo_info> algo_info_vec;
22 
23 std::map<std::string, JetAlgorithm> algorithms;
24 
25 struct loaderObj
26 {
27  loaderObj() {
28  static bool init = false;
29  if (!init)
30  {
31  algorithms["KT"] = kt_algorithm;
32  algorithms["CAMBRIDGE"] = cambridge_algorithm;
33  algorithms["ANTIKT"] = antikt_algorithm;
34  algorithms["GENKT"] = genkt_algorithm;
35  algorithms["CAMBRIDGE_FOR_PASSIVE"] = cambridge_for_passive_algorithm;
36  algorithms["GENKT_FOR_PASSIVE"] = genkt_for_passive_algorithm;
37  algorithms["EE_KT"] = ee_kt_algorithm;
38  algorithms["EE_GENKT"] = ee_genkt_algorithm;
39  algorithms["PLUGIN"] = plugin_algorithm ;
40  }
41  }
42 };
43 
44 loaderObj loader;
45 
46 void
47 hijfst_control(int enable, vector<string> valgorithm, vector<float> vR, vector<int> vPID, vector<float> vEtaMin, vector<float> vEtaMax, vector<float> vEtMin)
48 {
49  enablep = (enable==1) ? true: false;
50 
51  algo_info_vec.clear();
52  for (unsigned int i = 0; i < valgorithm.size(); ++i)
53  {
54  string algorithmName = valgorithm[i];
55  transform(algorithmName.begin(), algorithmName.end(), algorithmName.begin(), ::toupper);
56  algo_info algo;
57  algo.algorithm = ((algorithms.find(algorithmName) == algorithms.end()) ? antikt_algorithm : algorithms[algorithmName]);
58  algo.R = vR[i];
59  algo.PID = vPID[i];
60  algo.EtaMin = vEtaMin[i];
61  algo.EtaMax = vEtaMax[i];
62  algo.EtMin = vEtMin[i];
63  algo_info_vec.push_back(algo);
64  }
65 }
66 
67 extern "C"
68 void
69 hijfst_(int *n, int *N, int *K, float *P, float *V)
70 {
71  if (!enablep) return;
72  const int M = *N;
73 
74  int *K1 = new(K) int[M];
75  int *K2 = new(K+M) int[M];
76  int *K3 = new(K+2*M) int[M];
77  int *K4 = new(K+3*M) int[M];
78  int *K5 = new(K+4*M) int[M];
79 
80  float *px = new(P) float[M];
81  float *py = new(P+M) float[M];
82  float *pz = new(P+2*M) float[M];
83  float *E = new(P+3*M) float[M];
84  float *m = new(P+4*M) float[M];
85 
86  float *V1 = new(V) float[M];
87  float *V2 = new(V+M) float[M];
88  float *V3 = new(V+2*M) float[M];
89  float *V4 = new(V+3*M) float[M];
90  float *V5 = new(V+4*M) float[M];
91 
92  std::vector<PseudoJet> particles;
93 
94  for (int i = 0; i < *n; i++)
95  {
96  // We only want real, final state particles
97  if (K1[i] == 1)
98  {
99  particles.push_back(PseudoJet(px[i], py[i], pz[i], E[i]));
100  }
101  }
102 
103  for (vector<algo_info>::iterator it = algo_info_vec.begin(); it != algo_info_vec.end(); ++it)
104  {
105  algo_info a = *it;
106 
107  // Set the algorithm and R value
108  JetDefinition jet_def(a.algorithm, a.R);
109 
110  ClusterSequence cs(particles, jet_def);
111  std::vector<PseudoJet> jets = cs.inclusive_jets();
112 
113  // Loop over all the "jets" and add their kinematic properties to
114  // the LUJET common block. Set their status to 103 to distinguish
115  // them.
116  for (unsigned i = 0; i < jets.size(); i++)
117  {
118  // These cuts should be configurable!
119 // if (abs(jets[i].eta()) > 2.0 or jets[i].E() < 5.0)
120  if (jets[i].eta() < a.EtaMin or jets[i].eta() > a.EtaMax or jets[i].Et() < a.EtMin)
121  {
122  continue;
123  }
124  int j = (*n)++;
125 
126  K1[j] = 103;
127  K2[j] = a.PID;
128  K3[j] = 0;
129  K4[j] = 0;
130  K5[j] = 0;
131 
132  px[j] = jets[i].px();
133  py[j] = jets[i].py();
134  pz[j] = jets[i].pz();
135  E[j] = jets[i].E();
136  m[j] = jets[i].m();
137 
138  V1[j] = 0.0;
139  V2[j] = 0.0;
140  V3[j] = 0.0;
141  V4[j] = 0.0;
142  V5[j] = 0.0;
143  }
144  }
145 }
146