Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PJTranslator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PJTranslator.cc
1 //----------------------------------------------------------
2 // Module for flow jet reconstruction in sPHENIX
3 // J. Orjuela-Koop
4 // May 5 2015
5 //----------------------------------------------------------
6 
7 #include"PJTranslator.h"
8 #include"PseudoJetPlus.h"
9 //#include <fastjet/PseudoJet.hh>
10 
11 #include<calobase/RawCluster.h>
12 #include<calobase/RawClusterContainer.h>
13 #include <calobase/RawClusterUtility.h>
14 
15 #include <g4vertex/GlobalVertexMap.h>
16 #include <g4vertex/GlobalVertex.h>
17 
18 #include <g4hough/SvtxTrackMap.h>
19 #include <g4hough/SvtxTrack.h>
20 
21 #include<phool/PHCompositeNode.h>
22 #include<phool/PHNodeIterator.h>
23 #include<phool/PHNodeReset.h>
24 #include<fun4all/getClass.h>
25 
26 // #include <fastjet/JetDefinition.hh>
27 // #include <fastjet/ClusterSequence.hh>
28 // #include <fastjet/SISConePlugin.hh>
29 
30 #include <TF1.h>
31 // #include <TFile.h>
32 // #include <TH1.h>
33 // #include <TH2.h>
34 // #include <TLorentzVector.h>
35 // #include <TNtuple.h>
36 
37 // #include <cmath>
38 using namespace std;
39 
40 /*
41  * Constructor
42  */
43 PJTranslator::PJTranslator(const std::string &name):
44  SubsysReco(name)
45 {
46  //Define tolerance limits for track-cluster matching
47  match_tolerance_low = new TF1("match_tolerance_low","pol4");
48  match_tolerance_low->SetParameter(0,-0.470354);
49  match_tolerance_low->SetParameter(1,0.928888);
50  match_tolerance_low->SetParameter(2,-0.0958367);
51  match_tolerance_low->SetParameter(3,0.00748122);
52  match_tolerance_low->SetParameter(4,-0.000177858);
53 
54  match_tolerance_high = new TF1("match_tolerance_high","pol2");
55  match_tolerance_high->SetParameter(0,0.457184);
56  match_tolerance_high->SetParameter(1,1.24821);
57  match_tolerance_high->SetParameter(2,-0.00848157);
58 }
59 
60 /*
61  * Empty destructor
62  */
63 PJTranslator::~PJTranslator()
64 {
65 }
66 
67 /*
68  * Initialize module
69  */
70 int PJTranslator::Init(PHCompositeNode* topNode)
71 {
72  cout << "------PJTranslator::Init(PHCompositeNode*)------" << endl;
73  PHNodeIterator iter(topNode);
74 
76 
77  // // Looking for the DST node
78  // PHCompositeNode *dstNode
79  // = static_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode","DST"));
80  // if (!dstNode) {
81  // cout << PHWHERE << "DST Node missing, doing nothing." << endl;
82  // return Fun4AllReturnCodes::ABORTRUN;
83  // }
84 
85  particles = new PJContainer;
86 
87  // Create the PJ node if required
88  PHDataNode< PJContainer >* pjNode
89  = dynamic_cast< PHDataNode<PJContainer>* >(iter.findFirst("PHDataNode","PJ"));
90  if (!pjNode) {
91  pjNode = new PHDataNode< PJContainer >(particles, "PJ");
92  topNode->addNode(pjNode);
93  }
94 
95  return 0;
96 }
97 
98 /*
99  * Process event
100  */
101 int PJTranslator::process_event(PHCompositeNode* topNode)
102 {
103  cout << "------PJTranslator::process_event(PHCompositeNode*)------" << endl;
104 
105  if ( particles ) {
106  particles->data.clear();
107  } else {
108  cerr << " This should never happen." << endl;
109  return Fun4AllReturnCodes::ABORTRUN;
110 
111  }
112 
113  GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode,"GlobalVertexMap");
114  if (!vertexmap) {
115 
116  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;
117  assert(vertexmap); // force quit
118 
119  return Fun4AllReturnCodes::ABORTRUN;
120  }
121 
122  if (vertexmap->empty()) {
123  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;
124  return Fun4AllReturnCodes::ABORTRUN;
125  }
126 
127  GlobalVertex* vtx = vertexmap->begin()->second;
128  CLHEP::Hep3Vector vertex ;
129  if (vtx) vertex . set(vtx->get_x(),vtx->get_y(),vtx->get_z());
130 
131 
132  // vector<fastjet::PseudoJet> particles;
133  //-------------------------------------------------
134  // Gather constituents
135  //-------------------------------------------------
136 
137  //Loop over EMCAL clusters
138  //------------------------
139  RawClusterContainer *emc_clusters = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
140  for(unsigned int i=0; i<emc_clusters->size(); i++) {
141  RawCluster* part = emc_clusters->getCluster(i);
142 // double eta = part->get_eta();
143 // double phi = part->get_phi();
144 // double energy = part->get_energy()/sfEMCAL;
145 // double eT = energy/cosh(eta);
146 // double pz = eT*sinh(eta);
147 
148  CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*part, vertex);
149  E_vec_cluster = E_vec_cluster/sfEMCAL;
150  double eT = E_vec_cluster.perp();
151  double pz = E_vec_cluster.z();
152 
153 
154  if(eT<0.000001) {
155  eT = 0.001;
156  pz = eT*sinh(eta);
157  energy = sqrt(eT*eT + pz*pz);
158  }
161 
162  particles->data.push_back( fastjet::PseudoJet(eT*cos(phi),eT*sin(phi),pz,energy) );
163  }
164 
165 
166  //Loop over HCALIN clusters
167  //-------------------------
168  RawClusterContainer *hci_clusters = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_HCALIN");
169  for(unsigned int i=0; i<hci_clusters->size(); i++) {
170  RawCluster* part = hci_clusters->getCluster(i);
171 // double eta = part->get_eta();
172 // double phi = part->get_phi();
173 // double energy = part->get_energy()/sfHCALIN;
174 // double eT = energy/cosh(eta);
175 // double pz = eT*sinh(eta);
176 
177  CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*part, vertex);
178  E_vec_cluster = E_vec_cluster/sfHCALIN;
179  double eT = E_vec_cluster.perp();
180  double pz = E_vec_cluster.z();
181 
182  if(eT<0.000001) {
183  eT = 0.001;
184  pz = eT*sinh(eta);
185  energy = sqrt(eT*eT + pz*pz);
186  }
187 
190  particles->data.push_back( fastjet::PseudoJet(eT*cos(phi),eT*sin(phi),pz,energy) );
191  }
192 
193  //Loop over HCALOUT clusters
194  //-------------------------
195  RawClusterContainer *hco_clusters = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_HCALOUT");
196  for(unsigned int i=0; i<hco_clusters->size(); i++) {
197  RawCluster* part = hco_clusters->getCluster(i);
198 // double eta = part->get_eta();
199 // double phi = part->get_phi();
200 // double energy = part->get_energy()/sfHCALOUT;
201 // double eT = energy/cosh(eta);
202 // double pz = eT*sinh(eta);
203 
204  CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*part, vertex);
205  E_vec_cluster = E_vec_cluster/sfHCALOUT;
206  double eT = E_vec_cluster.perp();
207  double pz = E_vec_cluster.z();
208 
209  if(eT<0.000001) {
210  eT = 0.001;
211  pz = eT*sinh(eta);
212  energy = sqrt(eT*eT + pz*pz);
213  }
214 
217  particles->data.push_back( fastjet::PseudoJet(eT*cos(phi),eT*sin(phi),pz,energy) );
218  }
219 
220  //Get reconstructed tracks from nodetree
221  // -------------------------------------
222  // const float M = 0.139; //Assume pion mass
223  const float M = 0;
224 
225  SvtxTrackMap* reco_tracks = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
226  for(SvtxTrackMap::Iter iter = reco_tracks->begin(); iter != reco_tracks->end(); ++iter) {
227  SvtxTrack *trk = &iter->second;
228  double px = trk->get3Momentum(0);
229  double py = trk->get3Momentum(1);
230  double pz = trk->get3Momentum(2);
231  // double pt = sqrt(px*px + py*py);
232  double p = sqrt(px*px + py*py + pz*pz);
233  double track_energy = sqrt(p*p + M*M); //Assume pion mass
234  double phi = atan2(py,px);
235  double eta = -log(tan(acos(pz/p)/2.0));
236  double et = track_energy/cosh(eta);
237 
238  //Account for angle wrap-around
239  while( phi < 0 ) { phi += 2.0*M_PI; }
240  while( phi >= 2*M_PI ) { phi -= 2.0*M_PI; }
241 
242  //Quality cut on tracks
243  if(trk->getQuality() > 3.0) continue;
244 
245  // //Find ID of clusters that match to track in each layer
246  // int emcID = -1;
247  // int hciID = -1;
248  // int hcoID = -1;
249  // emcID = (int)trk->get_cal_cluster_id(SvtxTrack::CEMC);
250  // hciID = (int)trk->get_cal_cluster_id(SvtxTrack::HCALIN);
251  // hcoID = (int)trk->get_cal_cluster_id(SvtxTrack::HCALOUT);
252 
253  // //Find energy deposited by track in each layer
254  // tlvmap::iterator it = emc_map.find(emcID);
255  // if(it != emc_map.end()) {
256  // emc_energy = emc_map[emcID]->Energy();
257  // }
258 
259  // it = hci_map.find(hciID);
260  // if(it != hci_map.end()) {
261  // hci_energy = hci_map[hciID]->Energy();
262  // }
263 
264  // it = hco_map.find(hcoID);
265  // if(it != hco_map.end()) {
266  // hco_energy = hco_map[hcoID]->Energy();
267  // }
268 
269  // cluster_energy = emc_energy + hci_energy + hco_energy;
270 
271  // //Does the track match the cluster to within tolerance?
272  // // *matched = 0 --> clus_energy < track_energy
273  // // *matched = 1 --> clus_energy > track_energy
274  // // *matched = 2 --> clus_energy = track_energy
275  // int matched = -1;
276  // matched = get_matched(cluster_energy,track_energy);
277 
278  // //If matched = 1, remove track energy from clusters
279  // if(matched == 1) {
280  // float fracEnergyEMC = emc_energy/cluster_energy;
281  // float fracEnergyHCI = hci_energy/cluster_energy;
282  // float fracEnergyHCO = hco_energy/cluster_energy;
283 
284  // it = emc_map.find(emcID);
285  // if(it!=emc_map.end()) {
286  // (emc_map.find(emcID)->second)->SetE(emc_energy - fracEnergyEMC*track_energy);
287  // }
288 
289  // it = hci_map.find(hciID);
290  // if(it!=hci_map.end()) {
291  // (hci_map.find(hciID)->second)->SetE(hci_energy - fracEnergyHCI*track_energy);
292  // }
293 
294  // it = hco_map.find(hcoID);
295  // if(it!=hco_map.end()) {
296  // (hco_map.find(hcoID)->second)->SetE(hco_energy - fracEnergyHCO*track_energy);
297  // }
298  // } else if(matched == 2) {
299  // it = emc_map.find(emcID);
300  // if(it!=emc_map.end()) {
301  // delete emc_map[emcID];
302  // emc_map.erase(emcID);
303  // }
304 
305  // it = hci_map.find(hciID);
306  // if(it!=emc_map.end()) {
307  // delete hci_map[hciID];
308  // hci_map.erase(hciID);
309  // }
310 
311  // it = hco_map.find(hcoID);
312  // if(it!=hco_map.end()) {
313  // delete hco_map[hcoID];
314  // hco_map.erase(hcoID);
315  // }
316  // } else if(matched == 0) {
317  // continue;
318  // }
319 
320  // //Add perfectly matched and partially matched tracks to flow particle container
321  // if(et<0.000001) {
322  // et = 0.001;
323  // pz = et*sinh(eta);
324  // track_energy = sqrt(et*et + pz*pz);
325  // }
326 
327  particles->data.push_back( fastjet::PseudoJet (et*cos(phi),et*sin(phi),pz,track_energy) );
328  }
329 
330 
331  cout << "EMC: " << emc_clusters->size() << endl;
332  cout << "HCO: " << hco_clusters->size() << endl;
333  cout << "HCI: " << hci_clusters->size() << endl;
334  cout << "Tracker: " << reco_tracks->size() << endl;
335  cout << particles->data.size() << endl;
336 
337 
338  return 0;
339 }
340 
341 /*
342  * Given a track and cluster energies, determine if they match within tolerance
343  */
344 int PJTranslator::get_matched(double clus_energy, double track_energy)
345 {
346  double limLo = match_tolerance_low->Eval(track_energy);
347  double limHi = match_tolerance_high->Eval(track_energy);
348 
349  int matched = -1;
350 
351  if(clus_energy < limLo)
352  {
353  //Track energy greater than cluster energy. Most likely a fake
354  matched = 0;
355  }
356  else if(clus_energy > limHi)
357  {
358  //Contaminated cluster
359  matched = 1;
360  }
361  else
362  {
363  //Track and cluster match to within reason
364  matched = 2;
365  }
366 
367  return matched;
368 }
369 
370 /*
371  * Write jets to node tree
372  */
373 /*
374 todo - we are working on a jet object
375 void PHAJMaker::save_jets_to_nodetree()
376 {
377  flow_jet_container = new PHPyJetContainerV2();
378 }
379 */