sPhenixAj/src/PJTranslator.cc

0001 //----------------------------------------------------------
0002 // Module for flow jet reconstruction in sPHENIX
0003 // J. Orjuela-Koop
0004 // May 5 2015
0005 //----------------------------------------------------------
0006
0007 #include"PJTranslator.h"
0008 #include"PseudoJetPlus.h"
0009 //#include <fastjet/PseudoJet.hh>
0010
0011 #include<calobase/RawCluster.h>
0012 #include<calobase/RawClusterContainer.h>
0013 #include <calobase/RawClusterUtility.h>
0014
0015 #include <g4vertex/GlobalVertexMap.h>
0016 #include <g4vertex/GlobalVertex.h>
0017
0018 #include <g4hough/SvtxTrackMap.h>
0019 #include <g4hough/SvtxTrack.h>
0020
0021 #include<phool/PHCompositeNode.h>
0022 #include<phool/PHNodeIterator.h>
0023 #include<phool/PHNodeReset.h>
0024 #include<fun4all/getClass.h>
0025
0026 // #include <fastjet/JetDefinition.hh>
0027 // #include <fastjet/ClusterSequence.hh>
0028 // #include <fastjet/SISConePlugin.hh>
0029
0030 #include <TF1.h>
0031 // #include <TFile.h>
0032 // #include <TH1.h>
0033 // #include <TH2.h>
0034 // #include <TLorentzVector.h>
0035 // #include <TNtuple.h>
0036
0037 // #include <cmath>
0038 using namespace std;
0039
0040 /*
0041  * Constructor
0042  */
0043 PJTranslator::PJTranslator(const std::string &name):
0044   SubsysReco(name)
0045 {
0046   //Define tolerance limits for track-cluster matching
0047   match_tolerance_low = new TF1("match_tolerance_low","pol4");
0048   match_tolerance_low->SetParameter(0,-0.470354);
0049   match_tolerance_low->SetParameter(1,0.928888);
0050   match_tolerance_low->SetParameter(2,-0.0958367);
0051   match_tolerance_low->SetParameter(3,0.00748122);
0052   match_tolerance_low->SetParameter(4,-0.000177858);
0053
0054   match_tolerance_high = new TF1("match_tolerance_high","pol2");
0055   match_tolerance_high->SetParameter(0,0.457184);
0056   match_tolerance_high->SetParameter(1,1.24821);
0057   match_tolerance_high->SetParameter(2,-0.00848157);
0058 }
0059
0060 /*
0061  * Empty destructor
0062  */
0063 PJTranslator::~PJTranslator()
0064 {
0065 }
0066
0067 /*
0068  * Initialize module
0069  */
0070 int PJTranslator::Init(PHCompositeNode* topNode)
0071 {
0072   cout << "------PJTranslator::Init(PHCompositeNode*)------" << endl;
0073   PHNodeIterator iter(topNode);
0074
0075   /// \TODO: Create your own composite node, then attach to it.
0076
0077   // // Looking for the DST node
0078   // PHCompositeNode *dstNode
0079   //   = static_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode","DST"));
0080   // if (!dstNode) {
0081   //   cout << PHWHERE << "DST Node missing, doing nothing." << endl;
0082   //   return Fun4AllReturnCodes::ABORTRUN;
0083   // }
0084
0085   particles = new PJContainer;
0086
0087   // Create the PJ node if required
0088   PHDataNode< PJContainer >* pjNode
0089     = dynamic_cast< PHDataNode<PJContainer>* >(iter.findFirst("PHDataNode","PJ"));
0090   if (!pjNode) {
0091     pjNode = new PHDataNode< PJContainer >(particles, "PJ");
0092     topNode->addNode(pjNode);
0093   }
0094
0095   return 0;
0096 }
0097
0098 /*
0099  * Process event
0100  */
0101 int PJTranslator::process_event(PHCompositeNode* topNode)
0102 {
0103   cout << "------PJTranslator::process_event(PHCompositeNode*)------" << endl;
0104
0105   if ( particles ) {
0106     particles->data.clear();
0107   } else {
0108     cerr << " This should never happen." << endl;
0109     return Fun4AllReturnCodes::ABORTRUN;
0110
0111   }
0112
0113   GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode,"GlobalVertexMap");
0114   if (!vertexmap) {
0115
0116     cout <<"PJTranslator::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;
0117     assert(vertexmap); // force quit
0118
0119     return Fun4AllReturnCodes::ABORTRUN;
0120   }
0121
0122   if (vertexmap->empty()) {
0123     cout <<"PJTranslator::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;
0124     return Fun4AllReturnCodes::ABORTRUN;
0125   }
0126
0127   GlobalVertex* vtx = vertexmap->begin()->second;
0128   CLHEP::Hep3Vector vertex ;
0129   if (vtx) vertex . set(vtx->get_x(),vtx->get_y(),vtx->get_z());
0130
0131
0132   //  vector<fastjet::PseudoJet> particles;
0133   //-------------------------------------------------
0134   // Gather constituents
0135   //-------------------------------------------------
0136
0137   //Loop over EMCAL clusters
0138   //------------------------
0139  RawClusterContainer *emc_clusters = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
0140   for(unsigned int i=0; i<emc_clusters->size(); i++) {
0141     RawCluster* part = emc_clusters->getCluster(i);
0142 //    double eta = part->get_eta();
0143 //    double phi = part->get_phi();
0144 //    double energy = part->get_energy()/sfEMCAL;
0145 //    double eT = energy/cosh(eta);
0146 //    double pz = eT*sinh(eta);
0147
0148     CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*part, vertex);
0149     E_vec_cluster = E_vec_cluster/sfEMCAL;
0150     double eT = E_vec_cluster.perp();
0151     double pz = E_vec_cluster.z();
0152
0153
0154     if(eT<0.000001) {
0155     eT = 0.001;
0156     pz = eT*sinh(eta);
0157     energy = sqrt(eT*eT + pz*pz);
0158       }
0159     /// \TODO For EMCAL jets
0160     /// \TODO - Add kinematic cuts, add quality cuts, - add corrections
0161
0162     particles->data.push_back( fastjet::PseudoJet(eT*cos(phi),eT*sin(phi),pz,energy) );
0163   }
0164
0165
0166   //Loop over HCALIN clusters
0167   //-------------------------
0168   RawClusterContainer *hci_clusters = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_HCALIN");
0169   for(unsigned int i=0; i<hci_clusters->size(); i++) {
0170     RawCluster* part = hci_clusters->getCluster(i);
0171 //    double eta = part->get_eta();
0172 //    double phi = part->get_phi();
0173 //    double energy = part->get_energy()/sfHCALIN;
0174 //    double eT = energy/cosh(eta);
0175 //    double pz = eT*sinh(eta);
0176
0177     CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*part, vertex);
0178     E_vec_cluster = E_vec_cluster/sfHCALIN;
0179     double eT = E_vec_cluster.perp();
0180     double pz = E_vec_cluster.z();
0181
0182     if(eT<0.000001)  {
0183       eT = 0.001;
0184       pz = eT*sinh(eta);
0185       energy = sqrt(eT*eT + pz*pz);
0186     }
0187
0188     /// \TODO For HCALIN jets
0189     /// \TODO - Add kinematic cuts, add quality cuts, - add corrections
0190     particles->data.push_back( fastjet::PseudoJet(eT*cos(phi),eT*sin(phi),pz,energy) );
0191   }
0192
0193   //Loop over HCALOUT clusters
0194   //-------------------------
0195   RawClusterContainer *hco_clusters = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_HCALOUT");
0196   for(unsigned int i=0; i<hco_clusters->size(); i++) {
0197     RawCluster* part = hco_clusters->getCluster(i);
0198 //    double eta = part->get_eta();
0199 //    double phi = part->get_phi();
0200 //    double energy = part->get_energy()/sfHCALOUT;
0201 //    double eT = energy/cosh(eta);
0202 //    double pz = eT*sinh(eta);
0203
0204     CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*part, vertex);
0205     E_vec_cluster = E_vec_cluster/sfHCALOUT;
0206     double eT = E_vec_cluster.perp();
0207     double pz = E_vec_cluster.z();
0208
0209     if(eT<0.000001) {
0210       eT = 0.001;
0211       pz = eT*sinh(eta);
0212       energy = sqrt(eT*eT + pz*pz);
0213     }
0214
0215     /// \TODO For HCALOUT jets
0216     /// \TODO - Add kinematic cuts, add quality cuts, - add corrections
0217     particles->data.push_back( fastjet::PseudoJet(eT*cos(phi),eT*sin(phi),pz,energy) );
0218   }
0219
0220   //Get reconstructed tracks from nodetree
0221   // -------------------------------------
0222   //  const float M = 0.139; //Assume pion mass
0223   const float M = 0;
0224
0225   SvtxTrackMap* reco_tracks = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
0226   for(SvtxTrackMap::Iter iter = reco_tracks->begin(); iter != reco_tracks->end(); ++iter) {
0227     SvtxTrack *trk = &iter->second;
0228     double px = trk->get3Momentum(0);
0229     double py = trk->get3Momentum(1);
0230     double pz = trk->get3Momentum(2);
0231     // double pt = sqrt(px*px + py*py);
0232     double p = sqrt(px*px + py*py + pz*pz);
0233     double track_energy = sqrt(p*p + M*M); //Assume pion mass
0234     double phi = atan2(py,px);
0235     double eta = -log(tan(acos(pz/p)/2.0));
0236     double et = track_energy/cosh(eta);
0237
0238     //Account for angle wrap-around
0239     while( phi < 0 )       {  phi += 2.0*M_PI; }
0240     while( phi >= 2*M_PI ) {  phi -= 2.0*M_PI; }
0241
0242     //Quality cut on tracks
0243     if(trk->getQuality() > 3.0) continue;
0244
0245     // //Find ID of clusters that match to track in each layer
0246     // int emcID = -1;
0247     // int hciID = -1;
0248     // int hcoID = -1;
0249     // emcID = (int)trk->get_cal_cluster_id(SvtxTrack::CEMC);
0250     // hciID = (int)trk->get_cal_cluster_id(SvtxTrack::HCALIN);
0251     // hcoID = (int)trk->get_cal_cluster_id(SvtxTrack::HCALOUT);
0252
0253     // //Find energy deposited by track in each layer
0254     // tlvmap::iterator it = emc_map.find(emcID);
0255     // if(it != emc_map.end()) {
0256     //   emc_energy = emc_map[emcID]->Energy();
0257     // }
0258
0259     // it = hci_map.find(hciID);
0260     // if(it != hci_map.end()) {
0261     //   hci_energy = hci_map[hciID]->Energy();
0262     // }
0263
0264     // it = hco_map.find(hcoID);
0265     // if(it != hco_map.end()) {
0266     //   hco_energy = hco_map[hcoID]->Energy();
0267     // }
0268
0269     // cluster_energy = emc_energy + hci_energy + hco_energy;
0270
0271     // //Does the track match the cluster to within tolerance?
0272     // //  *matched = 0 --> clus_energy < track_energy
0273     // //  *matched = 1 --> clus_energy > track_energy
0274     // //  *matched = 2 --> clus_energy = track_energy
0275     // int matched = -1;
0276     // matched = get_matched(cluster_energy,track_energy);
0277
0278     // //If matched = 1, remove track energy from clusters
0279     // if(matched == 1) {
0280     //   float fracEnergyEMC = emc_energy/cluster_energy;
0281     //   float fracEnergyHCI = hci_energy/cluster_energy;
0282     //   float fracEnergyHCO = hco_energy/cluster_energy;
0283
0284     //   it = emc_map.find(emcID);
0285     //   if(it!=emc_map.end()) {
0286     //  (emc_map.find(emcID)->second)->SetE(emc_energy - fracEnergyEMC*track_energy);
0287     //   }
0288
0289     //   it = hci_map.find(hciID);
0290     //   if(it!=hci_map.end()) {
0291     //  (hci_map.find(hciID)->second)->SetE(hci_energy - fracEnergyHCI*track_energy);
0292     //   }
0293
0294     //   it = hco_map.find(hcoID);
0295     //   if(it!=hco_map.end()) {
0296     //  (hco_map.find(hcoID)->second)->SetE(hco_energy - fracEnergyHCO*track_energy);
0297     //   }
0298     // } else if(matched == 2) {
0299     //  it = emc_map.find(emcID);
0300     //  if(it!=emc_map.end()) {
0301     //    delete emc_map[emcID];
0302     //    emc_map.erase(emcID);
0303     //  }
0304
0305     //  it = hci_map.find(hciID);
0306     //  if(it!=emc_map.end()) {
0307     //    delete hci_map[hciID];
0308     //    hci_map.erase(hciID);
0309     //  }
0310
0311     //  it = hco_map.find(hcoID);
0312     //  if(it!=hco_map.end()) {
0313     //    delete hco_map[hcoID];
0314     //    hco_map.erase(hcoID);
0315     //  }
0316     // } else if(matched == 0) {
0317     //   continue;
0318     // }
0319
0320     // //Add perfectly matched and partially matched tracks to flow particle container
0321     // if(et<0.000001) {
0322     //   et = 0.001;
0323     //   pz = et*sinh(eta);
0324     //   track_energy = sqrt(et*et + pz*pz);
0325     // }
0326
0327     particles->data.push_back( fastjet::PseudoJet (et*cos(phi),et*sin(phi),pz,track_energy) );
0328   }
0329
0330
0331   cout << "EMC:     " << emc_clusters->size() << endl;
0332   cout << "HCO:     " << hco_clusters->size() << endl;
0333   cout << "HCI:     " << hci_clusters->size() << endl;
0334   cout << "Tracker: " << reco_tracks->size() << endl;
0335   cout << particles->data.size() << endl;
0336
0337
0338   return 0;
0339 }
0340
0341 /*
0342  * Given a track and cluster energies, determine if they match within tolerance
0343  */
0344 int PJTranslator::get_matched(double clus_energy, double track_energy)
0345 {
0346   double limLo = match_tolerance_low->Eval(track_energy);
0347   double limHi = match_tolerance_high->Eval(track_energy);
0348
0349   int matched = -1;
0350
0351   if(clus_energy < limLo)
0352     {
0353       //Track energy greater than cluster energy. Most likely a fake
0354       matched = 0;
0355     }
0356   else if(clus_energy > limHi)
0357     {
0358       //Contaminated cluster
0359       matched = 1;
0360     }
0361   else
0362     {
0363       //Track and cluster match to within reason
0364       matched = 2;
0365     }
0366
0367   return matched;
0368 }
0369
0370 /*
0371  * Write jets to node tree
0372  */
0373 /*
0374 todo - we are working on a jet object
0375 void PHAJMaker::save_jets_to_nodetree()
0376 {
0377   flow_jet_container = new PHPyJetContainerV2();
0378 }
0379 */