hijing/src/hijfst.cc

0001 #include <fastjet/ClusterSequence.hh>
0002
0003 #include <algorithm>
0004 #include <map>
0005 #include <string>
0006 #include <vector>
0007
0008 bool enablep = false;
0009 struct algo_info
0010 {
0011   fastjet::JetAlgorithm algorithm;
0012   double R;
0013   int PID;
0014   double EtaMin;
0015   double EtaMax;
0016   double EtMin;
0017 };
0018 std::vector<algo_info> algo_info_vec;
0019
0020 std::map<std::string, fastjet::JetAlgorithm> algorithms;
0021
0022 struct loaderObj
0023 {
0024   loaderObj()
0025   {
0026     static bool init = false;
0027     if (!init)
0028     {
0029       algorithms["KT"] = fastjet::kt_algorithm;
0030       algorithms["CAMBRIDGE"] = fastjet::cambridge_algorithm;
0031       algorithms["ANTIKT"] = fastjet::antikt_algorithm;
0032       algorithms["GENKT"] = fastjet::genkt_algorithm;
0033       algorithms["CAMBRIDGE_FOR_PASSIVE"] = fastjet::cambridge_for_passive_algorithm;
0034       algorithms["GENKT_FOR_PASSIVE"] = fastjet::genkt_for_passive_algorithm;
0035       algorithms["EE_KT"] = fastjet::ee_kt_algorithm;
0036       algorithms["EE_GENKT"] = fastjet::ee_genkt_algorithm;
0037       algorithms["PLUGIN"] = fastjet::plugin_algorithm;
0038     }
0039   }
0040 };
0041
0042 loaderObj loader;
0043
0044 void hijfst_control(int enable, const std::vector<std::string> &valgorithm, const std::vector<float> &vR, const std::vector<int> &vPID, const std::vector<float> &vEtaMin, const std::vector<float> &vEtaMax, const std::vector<float> &vEtMin)
0045 {
0046   enablep = (enable == 1) ? true : false;
0047
0048   algo_info_vec.clear();
0049   for (unsigned int i = 0; i < valgorithm.size(); ++i)
0050   {
0051     std::string algorithmName = valgorithm[i];
0052     std::transform(algorithmName.begin(), algorithmName.end(), algorithmName.begin(), ::toupper);
0053     algo_info algo{};
0054     algo.algorithm = ((!algorithms.contains(algorithmName)) ? fastjet::antikt_algorithm : algorithms[algorithmName]);
0055     algo.R = vR[i];
0056     algo.PID = vPID[i];
0057     algo.EtaMin = vEtaMin[i];
0058     algo.EtaMax = vEtaMax[i];
0059     algo.EtMin = vEtMin[i];
0060     algo_info_vec.push_back(algo);
0061   }
0062 }
0063
0064 extern "C" void hijfst_(int *n, int *N, int *K, float *P, float *V)  // NOLINT(readability-non-const-parameter)
0065 {
0066   if (!enablep)
0067   {
0068     return;
0069   }
0070   const int M = *N;
0071
0072   int *K1 = new (K) int[M];
0073   int *K2 = new (K + M) int[M];
0074   int *K3 = new (K + 2 * M) int[M];
0075   int *K4 = new (K + 3 * M) int[M];
0076   int *K5 = new (K + 4 * M) int[M];
0077
0078   float *px = new (P) float[M];
0079   float *py = new (P + M) float[M];
0080   float *pz = new (P + 2 * M) float[M];
0081   float *E = new (P + 3 * M) float[M];
0082   float *m = new (P + 4 * M) float[M];
0083
0084   float *V1 = new (V) float[M];
0085   float *V2 = new (V + M) float[M];
0086   float *V3 = new (V + 2 * M) float[M];
0087   float *V4 = new (V + 3 * M) float[M];
0088   float *V5 = new (V + 4 * M) float[M];
0089
0090   std::vector<fastjet::PseudoJet> particles;
0091
0092   for (int i = 0; i < *n; i++)
0093   {
0094     // We only want real, final state particles
0095     // cppcheck-suppress uninitdata
0096     if (K1[i] == 1)
0097     {
0098       // NOLINTNEXTLINE(hicpp-use-emplace,modernize-use-emplace)
0099       particles.push_back(fastjet::PseudoJet(px[i], py[i], pz[i], E[i]));
0100     }
0101   }
0102
0103   for (auto a : algo_info_vec)
0104   {
0105     // Set the algorithm and R value
0106     fastjet::JetDefinition jet_def(a.algorithm, a.R);
0107
0108     fastjet::ClusterSequence cs(particles, jet_def);
0109     std::vector<fastjet::PseudoJet> jets = cs.inclusive_jets();
0110
0111     // Loop over all the "jets" and add their kinematic properties to
0112     // the LUJET common block.  Set their status to 103 to distinguish
0113     // them.
0114     for (const auto &jet : jets)
0115     {
0116       // These cuts should be configurable!
0117       //        if (abs(jets[i].eta()) > 2.0 or jets[i].E() < 5.0)
0118       if (jet.eta() < a.EtaMin || jet.eta() > a.EtaMax || jet.Et() < a.EtMin)
0119       {
0120         continue;
0121       }
0122       int j = (*n)++;
0123
0124       K1[j] = 103;
0125       K2[j] = a.PID;
0126       K3[j] = 0;
0127       K4[j] = 0;
0128       K5[j] = 0;
0129
0130       px[j] = jet.px();
0131       py[j] = jet.py();
0132       pz[j] = jet.pz();
0133       E[j] = jet.E();
0134       m[j] = jet.m();
0135
0136       V1[j] = 0.0;
0137       V2[j] = 0.0;
0138       V3[j] = 0.0;
0139       V4[j] = 0.0;
0140       V5[j] = 0.0;
0141     }
0142   }
0143 }