analysis/PHFlowJetMaker/PHFlowJetMaker.C

0001 //----------------------------------------------------------
0002 // Module for flow jet reconstruction in sPHENIX
0003 // J. Orjuela-Koop
0004 // May 5 2015
0005 //----------------------------------------------------------
0006
0007 #include <PHFlowJetMaker.h>
0008 #include <calobase/RawCluster.h>
0009 #include <calobase/RawClusterContainer.h>
0010 #include <calobase/RawClusterUtility.h>
0011
0012 #include <g4vertex/GlobalVertex.h>
0013 #include <g4vertex/GlobalVertexMap.h>
0014
0015 #include <g4hough/SvtxTrack.h>
0016 #include <g4hough/SvtxTrackMap.h>
0017
0018 #include <phool/PHCompositeNode.h>
0019 #include <phool/PHNodeIterator.h>
0020 #include <phool/PHNodeReset.h>
0021 #include <phool/PHObject.h>
0022 #include <phool/getClass.h>
0023
0024 #include <fun4all/Fun4AllReturnCodes.h>
0025
0026 #include <fastjet/ClusterSequence.hh>
0027 #include <fastjet/JetDefinition.hh>
0028 #include <fastjet/PseudoJet.hh>
0029 #include <fastjet/SISConePlugin.hh>
0030
0031 #include <g4jets/Jet.h>
0032 #include <g4jets/JetMapV1.h>
0033 #include <g4jets/JetV1.h>
0034
0035 #include <TF1.h>
0036 #include <TFile.h>
0037 #include <TH1.h>
0038 #include <TH2.h>
0039 #include <TMath.h>
0040 #include <TNtuple.h>
0041
0042 #include <iomanip>
0043 #include <iostream>
0044 #include <sstream>
0045
0046 using namespace std;
0047 using namespace Fun4AllReturnCodes;
0048
0049 typedef std::map<int, TLorentzVector*> tlvmap;
0050
0051 // Jin - this thing is obsolete. Clusters now already have calibrated energies
0052 const float PHFlowJetMaker::sfEMCAL = 0.03;
0053 const float PHFlowJetMaker::sfHCALIN = 0.071;
0054 const float PHFlowJetMaker::sfHCALOUT = 0.04;
0055
0056 /*
0057  * Constructor
0058  */
0059 PHFlowJetMaker::PHFlowJetMaker(const std::string& name, const std::string algorithm, double r_param)
0060   : SubsysReco(name)
0061 {
0062   flow_jet_map = NULL;
0063   this->algorithm = algorithm;
0064   this->r_param = r_param;
0065   //outfile = outputfile;
0066
0067   //Define parameters for jet reconstruction
0068   min_jet_pT = 1.0;
0069   fastjet::Strategy strategy = fastjet::Best;
0070
0071   if (algorithm == "AntiKt")
0072     fJetAlgorithm = new fastjet::JetDefinition(fastjet::antikt_algorithm, r_param, fastjet::E_scheme, strategy);
0073
0074   if (algorithm == "Kt")
0075     fJetAlgorithm = new fastjet::JetDefinition(fastjet::kt_algorithm, r_param, fastjet::E_scheme, strategy);
0076
0077   //Define tolerance limits for track-cluster matching
0078   match_tolerance_low = new TF1("match_tolerance_low", "pol4");
0079   match_tolerance_low->SetParameter(0, -0.470354);
0080   match_tolerance_low->SetParameter(1, 0.928888);
0081   match_tolerance_low->SetParameter(2, -0.0958367);
0082   match_tolerance_low->SetParameter(3, 0.00748122);
0083   match_tolerance_low->SetParameter(4, -0.000177858);
0084
0085   match_tolerance_high = new TF1("match_tolerance_high", "pol2");
0086   match_tolerance_high->SetParameter(0, 0.457184);
0087   match_tolerance_high->SetParameter(1, 1.24821);
0088   match_tolerance_high->SetParameter(2, -0.00848157);
0089 }
0090
0091 /*
0092  * Empty destructor
0093  */
0094 PHFlowJetMaker::~PHFlowJetMaker()
0095 {
0096 }
0097
0098 /*
0099  * Initialize module
0100  */
0101 int PHFlowJetMaker::Init(PHCompositeNode* topNode)
0102 {
0103   cout << "------PHFlowJetMaker::Init(PHCompositeNode*)------" << endl;
0104   create_node_tree(topNode);
0105
0106   return EVENT_OK;
0107 }
0108
0109 /*
0110  * Create node tree for flow jets
0111  */
0112 int PHFlowJetMaker::create_node_tree(PHCompositeNode* topNode)
0113 {
0114   PHNodeIterator iter(topNode);
0115
0116   //Get the DST node
0117   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0118   if (!dstNode)
0119   {
0120     cout << "DST Node missing. Doing nothing." << endl;
0121     return ABORTEVENT;
0122   }
0123
0124   //Get the JET node
0125   PHCompositeNode* jetNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "JETS"));
0126   if (!jetNode)
0127   {
0128     if (algorithm == "AntiKt")
0129       jetNode = new PHCompositeNode("ANTIKT");
0130
0131     else if (algorithm == "Kt")
0132       jetNode = new PHCompositeNode("KT");
0133
0134     else
0135       jetNode = new PHCompositeNode(algorithm);
0136
0137     dstNode->addNode(jetNode);
0138   }
0139
0140   PHNodeIterator jiter(jetNode);
0141
0142   PHCompositeNode* flowJetNode = dynamic_cast<PHCompositeNode*>(jiter.findFirst("PHCompositeNode", "FLOW_JETS"));
0143   if (!flowJetNode)
0144   {
0145     flowJetNode = new PHCompositeNode("FLOW_JETS");
0146     jetNode->addNode(flowJetNode);
0147   }
0148
0149   //Add flow jet node to jet node
0150   flow_jet_map = new JetMapV1();
0151   string nodeName = "Flow";
0152
0153   stringstream snodeName;
0154   snodeName << algorithm << "_" << nodeName << "_r" << setfill('0') << setw(2) << int(r_param * 10);
0155   nodeName = snodeName.str();
0156
0157   PHIODataNode<PHObject>* PHFlowJetNode = new PHIODataNode<PHObject>(flow_jet_map, nodeName.c_str(), "PHObject");
0158
0159   if (!flowJetNode->addNode(PHFlowJetNode))
0160   {
0161     cout << "PHFlowJetMaker::create_node_tree() - Can't add node to tree!" << endl;
0162   }
0163
0164   return EVENT_OK;
0165 }
0166
0167 /*
0168  * Process event
0169  */
0170 int PHFlowJetMaker::process_event(PHCompositeNode* topNode)
0171 {
0172   //-------------------------------------------------
0173   // Get Information from Node Tree
0174   //-------------------------------------------------
0175
0176   //Get calorimeter clusters from node tree
0177   RawClusterContainer* emc_clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
0178   RawClusterContainer* hci_clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_HCALIN");
0179   RawClusterContainer* hco_clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_HCALOUT");
0180
0181   //Get reconstructed tracks from nodetree
0182   SvtxTrackMap* reco_tracks = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0183
0184   //Vertex for converting cluster position to momentum direction
0185   GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0186   if (!vertexmap)
0187   {
0188     cout << "PHFlowJetMaker::process_event - Fatal Error - GlobalVertexMap node is missing. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << endl;
0189     assert(vertexmap);  // force quit
0190
0191     return Fun4AllReturnCodes::ABORTRUN;
0192   }
0193
0194   if (vertexmap->empty())
0195   {
0196     cout << "PHFlowJetMaker::process_event - Fatal Error - GlobalVertexMap node is empty. Please turn on the do_bbc or tracking reco flags in the main macro in order to reconstruct the global vertex." << endl;
0197
0198     return Fun4AllReturnCodes::ABORTRUN;
0199   }
0200   GlobalVertex* vtx = vertexmap->begin()->second;
0201   assert(vtx);
0202
0203   //-------------------------------------------------
0204   // Create Jets
0205   //-------------------------------------------------
0206
0207   //  //Create jets from raw clusters
0208   //  vector<fastjet::PseudoJet> raw_cluster_particles;
0209   //  create_calo_pseudojets(raw_cluster_particles, emc_clusters, hci_clusters, hco_clusters, vtx);
0210   //  fastjet::ClusterSequence jet_finder_raw(raw_cluster_particles, *fJetAlgorithm);
0211   //  vector<fastjet::PseudoJet> raw_cluster_jets = jet_finder_raw.inclusive_jets(min_jet_pT);
0212
0213   //Apply particle flow jets algorithm and create jets from flow particles
0214   vector<fastjet::PseudoJet> flow_particles;
0215   run_particle_flow(flow_particles, emc_clusters, hci_clusters, hco_clusters, reco_tracks, vtx);
0216   fastjet::ClusterSequence jet_finder_flow(flow_particles, *fJetAlgorithm);
0217   vector<fastjet::PseudoJet> flow_jets = jet_finder_flow.inclusive_jets(min_jet_pT);
0218
0219   //Create JetMap from flow jets
0220   if (algorithm == "AntiKt")
0221     flow_jet_map->set_algo(Jet::ANTIKT);
0222
0223   if (algorithm == "Kt")
0224     flow_jet_map->set_algo(Jet::KT);
0225
0226   flow_jet_map->set_par(r_param);
0227   flow_jet_map->insert_src(Jet::TRACK);
0228   flow_jet_map->insert_src(Jet::CEMC_CLUSTER);
0229   flow_jet_map->insert_src(Jet::HCALIN_CLUSTER);
0230   flow_jet_map->insert_src(Jet::HCALOUT_CLUSTER);
0231
0232   for (unsigned int i = 0; i < flow_jets.size(); i++)
0233   {
0234     JetV1* j = new JetV1();
0235     float px = flow_jets[i].px();
0236     float py = flow_jets[i].py();
0237     float pz = flow_jets[i].pz();
0238     float energy = flow_jets[i].E();
0239
0240     //      cout << "--px = " << px << endl;
0241
0242     j->set_px(px);
0243     j->set_py(py);
0244     j->set_pz(pz);
0245     j->set_e(energy);
0246
0247     flow_jet_map->insert(j);
0248   }
0249
0250   //  cout << "IDENTIFYING JET MAP" << endl;
0251   //  flow_jet_map->identify();
0252   //
0253   //  cout << "IDENTIFYING INDIVIDUAL JETS" << endl;
0254   //  for(JetMap::Iter it = flow_jet_map->begin(); it != flow_jet_map->end(); it++)
0255   //    {
0256   //      (it->second)->identify();
0257   //      cout << endl;
0258   //    }
0259   //
0260   //  cout << "FOUND " << flow_jet_map->size() << " FLOW JETS" << endl;
0261   return EVENT_OK;
0262 }
0263
0264 /*
0265  * Run particle flow algorithm
0266  */
0267 void PHFlowJetMaker::run_particle_flow(std::vector<fastjet::PseudoJet>& flow_particles, RawClusterContainer* emc_clusters, RawClusterContainer* hci_clusters, RawClusterContainer* hco_clusters, SvtxTrackMap* reco_tracks, GlobalVertex* vtx)
0268 {
0269   CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0270
0271   double emc_energy = 0;
0272   double hci_energy = 0;
0273   double hco_energy = 0;
0274   double cluster_energy = 0;
0275
0276   double px = 0;
0277   double py = 0;
0278   double pz = 0;
0279   //  double pt = 0;
0280   double et = 0;
0281   double p = 0;
0282   double track_energy = 0;
0283   double phi = 0;
0284   double eta = 0;
0285
0286   //Make track maps
0287   tlvmap emc_map;
0288   tlvmap hci_map;
0289   tlvmap hco_map;
0290
0291   for (unsigned int i = 0; i < emc_clusters->size(); i++)
0292   {
0293     RawCluster* part = emc_clusters->getCluster(i);
0294     //    double pT = (part->get_energy()) / cosh(part->get_eta());
0295
0296     CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*part, vertex);
0297
0298     emc_map[i] = new TLorentzVector();
0299     //    emc_map[i]->SetPtEtaPhiE(pT, part->get_eta(), part->get_phi(), part->get_energy());
0300     emc_map[i]->SetPxPyPzE(
0301         E_vec_cluster.x(),
0302         E_vec_cluster.y(),
0303         E_vec_cluster.z(),
0304         E_vec_cluster.mag());
0305   }
0306
0307   for (unsigned int i = 0; i < hci_clusters->size(); i++)
0308   {
0309     RawCluster* part = hci_clusters->getCluster(i);
0310     //    double pT = (part->get_energy()) / cosh(part->get_eta());
0311     CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*part, vertex);
0312     hci_map[i] = new TLorentzVector();
0313     hci_map[i]->
0314         //    SetPtEtaPhiE(pT, part->get_eta(), part->get_phi(), part->get_energy());
0315         SetPxPyPzE(
0316             E_vec_cluster.x(),
0317             E_vec_cluster.y(),
0318             E_vec_cluster.z(),
0319             E_vec_cluster.mag());
0320   }
0321
0322   for (unsigned int i = 0; i < hco_clusters->size(); i++)
0323   {
0324     RawCluster* part = hco_clusters->getCluster(i);
0325     //    double pT = (part->get_energy()) / cosh(part->get_eta());
0326     CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*part, vertex);
0327     hco_map[i] = new TLorentzVector();
0328     hco_map[i]->
0329         //    SetPtEtaPhiE(pT, part->get_eta(), part->get_phi(), part->get_energy());
0330         SetPxPyPzE(
0331             E_vec_cluster.x(),
0332             E_vec_cluster.y(),
0333             E_vec_cluster.z(),
0334             E_vec_cluster.mag());
0335   }
0336
0337   //Loop over all tracks
0338   for (SvtxTrackMap::Iter iter = reco_tracks->begin(); iter != reco_tracks->end(); ++iter)
0339   {
0340     SvtxTrack* trk = iter->second;
0341     px = trk->get_px();
0342     py = trk->get_py();
0343     pz = trk->get_pz();
0344     //      pt = sqrt(px*px + py*py);
0345     p = sqrt(px * px + py * py + pz * pz);
0346     track_energy = TMath::Sqrt(p * p + 0.139 * 0.139);  //Assume pion mass
0347     phi = atan2(py, px);
0348     eta = -log(tan(acos(pz / p) / 2.0));
0349     et = track_energy / cosh(eta);
0350
0351     //Account for angle wrap-around
0352     if (phi < 0)
0353     {
0354       phi = phi + 2 * TMath::Pi();
0355     }
0356     else if (phi > 2 * TMath::Pi())
0357     {
0358       phi = phi + 2 * TMath::Pi();
0359     }
0360
0361     //Quality cut on tracks
0362     if (trk->get_quality() > 3.0) continue;
0363
0364     //Find ID of clusters that match to track in each layer
0365     int emcID = -1;
0366     int hciID = -1;
0367     int hcoID = -1;
0368     emcID = (int) trk->get_cal_cluster_id(SvtxTrack::CEMC);
0369     hciID = (int) trk->get_cal_cluster_id(SvtxTrack::HCALIN);
0370     hcoID = (int) trk->get_cal_cluster_id(SvtxTrack::HCALOUT);
0371
0372     //Find energy deposited by track in each layer
0373     tlvmap::iterator it = emc_map.find(emcID);
0374     if (it != emc_map.end())
0375     {
0376       emc_energy = emc_map[emcID]->Energy();
0377     }
0378
0379     it = hci_map.find(hciID);
0380     if (it != hci_map.end())
0381     {
0382       hci_energy = hci_map[hciID]->Energy();
0383     }
0384
0385     it = hco_map.find(hcoID);
0386     if (it != hco_map.end())
0387     {
0388       hco_energy = hco_map[hcoID]->Energy();
0389     }
0390
0391     cluster_energy = emc_energy + hci_energy + hco_energy;
0392
0393     //Does the track match the cluster to within tolerance?
0394     //  *matched = 0 --> clus_energy < track_energy
0395     //  *matched = 1 --> clus_energy > track_energy
0396     //  *matched = 2 --> clus_energy = track_energy
0397     int matched = -1;
0398     matched = get_matched(cluster_energy, track_energy);
0399
0400     //If matched = 1, remove track energy from clusters
0401     if (matched == 1)
0402     {
0403       float fracEnergyEMC = emc_energy / cluster_energy;
0404       float fracEnergyHCI = hci_energy / cluster_energy;
0405       float fracEnergyHCO = hco_energy / cluster_energy;
0406
0407       it = emc_map.find(emcID);
0408       if (it != emc_map.end())
0409       {
0410         (emc_map.find(emcID)->second)->SetE(emc_energy - fracEnergyEMC * track_energy);
0411       }
0412
0413       it = hci_map.find(hciID);
0414       if (it != hci_map.end())
0415       {
0416         (hci_map.find(hciID)->second)->SetE(hci_energy - fracEnergyHCI * track_energy);
0417       }
0418
0419       it = hco_map.find(hcoID);
0420       if (it != hco_map.end())
0421       {
0422         (hco_map.find(hcoID)->second)->SetE(hco_energy - fracEnergyHCO * track_energy);
0423       }
0424     }
0425     else if (matched == 2)
0426     {
0427       it = emc_map.find(emcID);
0428       if (it != emc_map.end())
0429       {
0430         delete emc_map[emcID];
0431         emc_map.erase(emcID);
0432       }
0433
0434       it = hci_map.find(hciID);
0435       if (it != emc_map.end())
0436       {
0437         delete hci_map[hciID];
0438         hci_map.erase(hciID);
0439       }
0440
0441       it = hco_map.find(hcoID);
0442       if (it != hco_map.end())
0443       {
0444         delete hco_map[hcoID];
0445         hco_map.erase(hcoID);
0446       }
0447     }
0448     else if (matched == 0)
0449     {
0450       continue;
0451     }
0452
0453     //Add perfectly matched and partially matched tracks to flow particle container
0454     if (et < 0.000001)
0455     {
0456       et = 0.001;
0457       pz = et * sinh(eta);
0458       track_energy = sqrt(et * et + pz * pz);
0459     }
0460     fastjet::PseudoJet pseudoJet_track(et * cos(phi), et * sin(phi), pz, track_energy);
0461     flow_particles.push_back(pseudoJet_track);
0462   }
0463
0464   //Add remaining clusters to flow particle container
0465   for (tlvmap::iterator it = emc_map.begin(); it != emc_map.end(); it++)
0466   {
0467     double energy_clus = (it->second)->Energy();
0468     double eta_clus = (it->second)->Eta();
0469     double phi_clus = (it->second)->Phi();
0470     double et_clus = energy_clus / cosh(eta_clus);
0471     double pz_clus = et * sinh(eta_clus);
0472
0473     if (et_clus < 0.000001)
0474     {
0475       et_clus = 0.001;
0476       pz_clus = et_clus * sinh(eta_clus);
0477       energy_clus = sqrt(et_clus * et_clus + pz_clus * pz_clus);
0478     }
0479     fastjet::PseudoJet pseudoJet_clus(et_clus * cos(phi_clus), et_clus * sin(phi_clus), pz_clus, energy_clus);
0480     flow_particles.push_back(pseudoJet_clus);
0481   }
0482
0483   for (tlvmap::iterator it = hci_map.begin(); it != hci_map.end(); it++)
0484   {
0485     double energy_clus = (it->second)->Energy();
0486     double eta_clus = (it->second)->Eta();
0487     double phi_clus = (it->second)->Phi();
0488     double et_clus = energy_clus / cosh(eta_clus);
0489     double pz_clus = et * sinh(eta_clus);
0490
0491     if (et_clus < 0.000001)
0492     {
0493       et_clus = 0.001;
0494       pz_clus = et_clus * sinh(eta_clus);
0495       energy_clus = sqrt(et_clus * et_clus + pz_clus * pz_clus);
0496     }
0497     fastjet::PseudoJet pseudoJet_clus(et_clus * cos(phi_clus), et_clus * sin(phi_clus), pz_clus, energy_clus);
0498     flow_particles.push_back(pseudoJet_clus);
0499   }
0500
0501   for (tlvmap::iterator it = hco_map.begin(); it != hco_map.end(); it++)
0502   {
0503     double energy_clus = (it->second)->Energy();
0504     double eta_clus = (it->second)->Eta();
0505     double phi_clus = (it->second)->Phi();
0506     double et_clus = energy_clus / cosh(eta_clus);
0507     double pz_clus = et * sinh(eta_clus);
0508
0509     if (et_clus < 0.000001)
0510     {
0511       et_clus = 0.001;
0512       pz_clus = et_clus * sinh(eta_clus);
0513       energy_clus = sqrt(et_clus * et_clus + pz_clus * pz_clus);
0514     }
0515     fastjet::PseudoJet pseudoJet_clus(et_clus * cos(phi_clus), et_clus * sin(phi_clus), pz_clus, energy_clus);
0516     flow_particles.push_back(pseudoJet_clus);
0517   }
0518 }
0519
0520 /*
0521  * Create pseudojets from calorimeter clusters before applying particle flow algorithm
0522  * Jin - disable. Please use the standard cluster jet on the node tree
0523  */
0524 //void PHFlowJetMaker::create_calo_pseudojets(std::vector<fastjet::PseudoJet>& particles, RawClusterContainer* emc_clusters, RawClusterContainer* hci_clusters, RawClusterContainer* hco_clusters, GlobalVertex* vtx)
0525 //{
0526 //  CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0527 //
0528 //  //Loop over EMCAL clusters
0529 //  for (unsigned int i = 0; i < emc_clusters->size(); i++)
0530 //  {
0531 //    RawCluster* part = emc_clusters->getCluster(i);
0532 //    double eta = part->get_eta();
0533 //    double phi = part->get_phi();
0534 //    double energy = part->get_energy() / sfEMCAL;
0535 //    double eT = energy / cosh(eta);
0536 //    double pz = eT * sinh(eta);
0537 //
0538 //    if (eT < 0.000001)
0539 //    {
0540 //      eT = 0.001;
0541 //      pz = eT * sinh(eta);
0542 //      energy = sqrt(eT * eT + pz * pz);
0543 //    }
0544 //
0545 //    fastjet::PseudoJet pseudoJet(eT * cos(phi), eT * sin(phi), pz, energy);
0546 //    particles.push_back(pseudoJet);
0547 //  }
0548 //
0549 //  //Loop over HCALIN clusters
0550 //  for (unsigned int i = 0; i < hci_clusters->size(); i++)
0551 //  {
0552 //    RawCluster* part = hci_clusters->getCluster(i);
0553 //    double eta = part->get_eta();
0554 //    double phi = part->get_phi();
0555 //    double energy = part->get_energy() / sfHCALIN;
0556 //    double eT = energy / cosh(eta);
0557 //    double pz = eT * sinh(eta);
0558 //
0559 //    if (eT < 0.000001)
0560 //    {
0561 //      eT = 0.001;
0562 //      pz = eT * sinh(eta);
0563 //      energy = sqrt(eT * eT + pz * pz);
0564 //    }
0565 //
0566 //    fastjet::PseudoJet pseudoJet(eT * cos(phi), eT * sin(phi), pz, energy);
0567 //    particles.push_back(pseudoJet);
0568 //  }
0569 //
0570 //  //Loop over HCALOUT clusters
0571 //  for (unsigned int i = 0; i < hco_clusters->size(); i++)
0572 //  {
0573 //    RawCluster* part = hco_clusters->getCluster(i);
0574 //    double eta = part->get_eta();
0575 //    double phi = part->get_phi();
0576 //    double energy = part->get_energy() / sfHCALOUT;
0577 //    double eT = energy / cosh(eta);
0578 //    double pz = eT * sinh(eta);
0579 //
0580 //    if (eT < 0.000001)
0581 //    {
0582 //      eT = 0.001;
0583 //      pz = eT * sinh(eta);
0584 //      energy = sqrt(eT * eT + pz * pz);
0585 //    }
0586 //
0587 //    fastjet::PseudoJet pseudoJet(eT * cos(phi), eT * sin(phi), pz, energy);
0588 //    particles.push_back(pseudoJet);
0589 //  }
0590 //}
0591
0592 /*
0593  * Given a track and cluster energies, determine if they match within tolerance
0594  */
0595 int PHFlowJetMaker::get_matched(double clus_energy, double track_energy)
0596 {
0597   double limLo = match_tolerance_low->Eval(track_energy);
0598   double limHi = match_tolerance_high->Eval(track_energy);
0599
0600   int matched = -1;
0601
0602   if (clus_energy < limLo)
0603   {
0604     //Track energy greater than cluster energy. Most likely a fake
0605     matched = 0;
0606   }
0607   else if (clus_energy > limHi)
0608   {
0609     //Contaminated cluster
0610     matched = 1;
0611   }
0612   else
0613   {
0614     //Track and cluster match to within reason
0615     matched = 2;
0616   }
0617
0618   return matched;
0619 }
0620
0621 /*
0622  * End
0623  */
0624 int PHFlowJetMaker::End(PHCompositeNode* topNode)
0625 {
0626   return EVENT_OK;
0627 }