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     static bool init = false;
0026     if (!init)
0027       {
0028     algorithms["KT"] = fastjet::kt_algorithm;
0029     algorithms["CAMBRIDGE"] = fastjet::cambridge_algorithm;
0030     algorithms["ANTIKT"] = fastjet::antikt_algorithm;
0031     algorithms["GENKT"] = fastjet::genkt_algorithm;
0032     algorithms["CAMBRIDGE_FOR_PASSIVE"] = fastjet::cambridge_for_passive_algorithm;
0033     algorithms["GENKT_FOR_PASSIVE"] = fastjet::genkt_for_passive_algorithm;
0034     algorithms["EE_KT"] = fastjet::ee_kt_algorithm;
0035     algorithms["EE_GENKT"] = fastjet::ee_genkt_algorithm;
0036     algorithms["PLUGIN"] = fastjet::plugin_algorithm ;
0037       }
0038   }
0039 };
0040
0041 loaderObj loader;
0042
0043 void
0044 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.find(algorithmName) == algorithms.end()) ? 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"
0065 void
0066 hijfst_(int *n, int *N, int *K, float *P, float *V)
0067 {
0068   if (!enablep) return;
0069   const int M = *N;
0070
0071   int *K1 = new(K)     int[M];
0072   int *K2 = new(K+M)   int[M];
0073   int *K3 = new(K+2*M) int[M];
0074   int *K4 = new(K+3*M) int[M];
0075   int *K5 = new(K+4*M) int[M];
0076
0077   float *px = new(P)     float[M];
0078   float *py = new(P+M)   float[M];
0079   float *pz = new(P+2*M) float[M];
0080   float *E  = new(P+3*M) float[M];
0081   float *m  = new(P+4*M) float[M];
0082
0083   float *V1 = new(V)     float[M];
0084   float *V2 = new(V+M)   float[M];
0085   float *V3 = new(V+2*M) float[M];
0086   float *V4 = new(V+3*M) float[M];
0087   float *V5 = new(V+4*M) float[M];
0088
0089   std::vector<fastjet::PseudoJet> particles;
0090
0091   for (int i = 0; i < *n; i++)
0092     {
0093       // We only want real, final state particles
0094  // cppcheck-suppress uninitdata
0095       if (K1[i] == 1)
0096     {
0097       particles.push_back(fastjet::PseudoJet(px[i], py[i], pz[i], E[i]));
0098     }
0099     }
0100
0101   for (std::vector<algo_info>::iterator it = algo_info_vec.begin(); it != algo_info_vec.end(); ++it)
0102     {
0103       algo_info a = *it;
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 (unsigned i = 0; i < jets.size(); i++)
0115     {
0116       // These cuts should be configurable!
0117 //        if (abs(jets[i].eta()) > 2.0 or jets[i].E() < 5.0)
0118       if (jets[i].eta() < a.EtaMin or jets[i].eta() > a.EtaMax or  jets[i].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] = jets[i].px();
0131       py[j] = jets[i].py();
0132       pz[j] = jets[i].pz();
0133       E[j]  = jets[i].E();
0134       m[j]  = jets[i].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 }
0144