Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FilterEventsUpsilon.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FilterEventsUpsilon.cc
1 
2 #include "FilterEventsUpsilon.h"
3 
4 #include <cstdlib>
5 #include <cstdio>
6 #include <iostream>
7 #include <iomanip>
8 #include <fstream>
9 
10 #include <TFile.h>
11 
12 #include <trackbase_historic/SvtxTrack.h>
13 #include <trackbase_historic/SvtxTrackMap.h>
14 #include <trackbase_historic/SvtxTrackMap_v2.h>
15 #include <trackbase_historic/SvtxVertex.h>
16 #include <trackbase_historic/SvtxVertexMap.h>
17 #include <trackbase_historic/TrackSeed.h>
18 #include <trackbase_historic/TrackSeed_v1.h>
19 #include <trackbase_historic/SvtxTrackSeed_v1.h>
20 #include <trackbase_historic/TrackSeedContainer.h>
21 #include <trackbase/TrkrDefs.h>
22 #include <trackbase/TrkrCluster.h>
23 #include <trackbase/TrkrClusterv4.h>
24 #include <trackbase/TrkrClusterContainer.h>
25 #include <g4vertex/GlobalVertexMap.h>
26 #include <g4vertex/GlobalVertex.h>
27 #include <calobase/RawClusterContainer.h>
28 #include <calobase/RawCluster.h>
29 #include <calobase/RawClusterv1.h>
30 
31 #include <phool/getClass.h>
32 #include <phool/recoConsts.h>
33 #include <phool/PHCompositeNode.h>
34 #include <phool/PHIODataNode.h>
35 #include <phool/PHNodeIterator.h>
36 #include <phool/PHRandomSeed.h>
37 #include <fun4all/Fun4AllReturnCodes.h>
38 
39 typedef PHIODataNode<PHObject> PHObjectNode_t;
40 
41 using namespace std;
42 
43 //==============================================================
44 
45 FilterEventsUpsilon::FilterEventsUpsilon(const std::string &name) : SubsysReco(name)
46 {
47  outnodename_trackmap = "SvtxTrackMap_ee";
48  outnodename_cemc_clusters = "CLUSTER_CEMC_ee";
49  EventNumber=0;
50  goodEventNumber=0;
51 }
52 
53 //-------------------------------------------------------------------------------
54 
55 int FilterEventsUpsilon::Init(PHCompositeNode *topNode)
56 {
57 
58 cout << "FilterEventsUpsilon::Init started..." << endl;
59 
60  PHNodeIterator iter(topNode);
61  //PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
62  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
63  if (!dstNode) { cerr << PHWHERE << " ERROR: DST node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
64 
65 // PHNodeIterator dstiter(dstNode);
66 // PHCompositeNode *svtxNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "SVTX"));
67 // if (!svtxNode) { cerr << PHWHERE << " ERROR: SVTX node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
68 //
69 // PHCompositeNode *cemcNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "CEMC"));
70 // if (!cemcNode) { cerr << PHWHERE << " ERROR: CEMC node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
71 
72  SvtxTrackMap_v2* trackmap = new SvtxTrackMap_v2();
73  PHCompositeNode *trackmapNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", outnodename_trackmap));
74  if (!trackmapNode)
75  {
76  PHObjectNode_t *trackmapNode = new PHIODataNode<PHObject>(trackmap,outnodename_trackmap.c_str(),"PHObject");
77  dstNode->addNode(trackmapNode);
78  //svtxNode->addNode(trackmapNode);
79  cout << PHWHERE << " INFO: added " << outnodename_trackmap << " node." << endl;
80  }
81  else { cout << PHWHERE << " INFO: " << outnodename_trackmap << " node already exists." << endl; }
82 
83  _cemc_clusters_ee = new RawClusterContainer();
84  PHIODataNode<PHObject> *clusterNode = new PHIODataNode<PHObject>(_cemc_clusters_ee, outnodename_cemc_clusters, "PHObject");
85  dstNode->addNode(clusterNode);
86  //cemcNode->addNode(clusterNode);
87 
88 // PHCompositeNode *cemcclusNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", outnodename_cemc_clusters));
89 // if (!cemcclusNode)
90 // {
91 // PHObjectNode_t *cemcclusNode = new PHIODataNode<PHObject>(cemc_clusters,outnodename_cemc_clusters.c_str(),"PHObject");
92 // dstNode->addNode(cemcclusNode);
93 // cout << PHWHERE << " INFO: added " << outnodename_cemc_clusters << " node." << endl;
94 // }
95 // else { cout << PHWHERE << " INFO: " << outnodename_cemc_clusters << " node already exists." << endl; }
96 
97 cout << "FilterEventsUpsilon::Init ended." << endl;
98  return Fun4AllReturnCodes::EVENT_OK;
99 }
100 
101 //--------------------------------------------------------------------------------
102 
103 int FilterEventsUpsilon::InitRun(PHCompositeNode *topNode)
104 {
105  return Fun4AllReturnCodes::EVENT_OK;
106 }
107 
108 //--------------------------------------------------------------------------------
109 
110 int FilterEventsUpsilon::GetNodes(PHCompositeNode *topNode)
111 {
112 
113  _topNode = topNode;
114 
115  //topNode->print();
116 
117  _trackmap_ee = findNode::getClass<SvtxTrackMap>(topNode, outnodename_trackmap);
118  if(!_trackmap_ee) { cerr << PHWHERE << "ERROR: Output SvtxTrackMap_ee node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
119 // else {cout << "Found SvtxTrackMap_ee node." << endl; }
120 
121  _cemc_clusters_ee = findNode::getClass<RawClusterContainer>(topNode, outnodename_cemc_clusters);
122  if(!_cemc_clusters_ee) { cerr << PHWHERE << "ERROR: CLUSTER_CEMC_ee node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
123 // else {cout << "Found CLUSTER_CEMC_ee node." << endl; }
124 
125  _trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
126  if(!_trackmap) { cerr << PHWHERE << "ERROR: SvtxTrackMap node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
127 
128 // _vtxmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
129 // if(!_vtxmap) { cout << "ERROR: SvtxVertexMap node not found!" << endl; return Fun4AllReturnCodes::ABORTEVENT; }
130 
131  _trackseedcontainer_svtx = findNode::getClass<TrackSeedContainer>(topNode, "SvtxTrackSeedContainer");
132  if(!_trackseedcontainer_svtx) { cerr << PHWHERE << "ERROR: SvtxTrackSeedContainer node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
133 
134  _trackseedcontainer_silicon = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
135  if(!_trackseedcontainer_silicon) { cerr << PHWHERE << "ERROR: SiliconTrackSeedContainer node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
136 
137  _trackseedcontainer_tpc = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
138  if(!_trackseedcontainer_tpc) { cerr << PHWHERE << "ERROR: TpcTrackSeedContainer node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
139 
140  _trkrclusters = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
141  if(!_trkrclusters) { cerr << PHWHERE << "ERROR: TRKR_CLUSTER node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
142 
143  _cemc_clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
144  if(!_cemc_clusters) { cerr << PHWHERE << "ERROR: CLUSTER_CEMC node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
145 
146  return Fun4AllReturnCodes::EVENT_OK;
147 }
148 
149 //--------------------------------------------------------------------
150 
151 TVector3 FilterEventsUpsilon::GetProjectionCEMC(SvtxTrack* track) {
152 
153  TVector3 projection; // 0,0,0
154 
155  vector<double> proj;
156  for (SvtxTrack::StateIter stateiter = track->begin_states(); stateiter != track->end_states(); ++stateiter)
157  {
158  SvtxTrackState *trackstate = stateiter->second;
159  if(trackstate) { proj.push_back(trackstate->get_pathlength()); }
160  }
161  double pathlength = proj[proj.size()-3]; // CEMC is next next to last, usually 93.5
162 
163  SvtxTrackState* trackstate = track->get_state(pathlength); // at CEMC inner face
164  if(trackstate) {
165  double track_x = trackstate->get_x();
166  double track_y = trackstate->get_y();
167  double track_z = trackstate->get_z();
168  projection.SetX(track_x);
169  projection.SetY(track_y);
170  projection.SetZ(track_z);
171  }
172 
173  return projection;
174 }
175 
176 //--------------------------------------------------------------------
177 
178 RawCluster* FilterEventsUpsilon::MatchClusterCEMC(SvtxTrack* track, RawClusterContainer* cemc_clusters, double &dphi, double &deta, double Zvtx) {
179 
180  RawCluster* returnCluster = NULL;
181  double track_eta = 99999.;
182  double track_phi = 99999.;
183  dphi = 99999.;
184  deta = 99999.;
185 
186  vector<double> proj;
187  for (SvtxTrack::StateIter stateiter = track->begin_states(); stateiter != track->end_states(); ++stateiter)
188  {
189  SvtxTrackState *trackstate = stateiter->second;
190  if(trackstate) { proj.push_back(trackstate->get_pathlength()); }
191  }
192  double pathlength = proj[proj.size()-3]; // CEMC is next next to last
193 
194  SvtxTrackState* trackstate = track->get_state(pathlength); // at CEMC inner face
195  if(trackstate) {
196  double track_x = trackstate->get_x();
197  double track_y = trackstate->get_y();
198  double track_z = trackstate->get_z() - Zvtx;
199  double track_r = sqrt(track_x*track_x+track_y*track_y);
200  track_eta = asinh( track_z / track_r );
201  track_phi = atan2( track_y, track_x );
202  } else { return returnCluster; }
203 
204  if(track_eta == 99999. || track_phi == 99999.) { return returnCluster; }
205  double dist = 99999.;
206 
207  RawClusterContainer::Range begin_end = cemc_clusters->getClusters();
208  RawClusterContainer::Iterator iter;
209  for (iter = begin_end.first; iter != begin_end.second; ++iter)
210  {
211  RawCluster* cluster = iter->second;
212  if(!cluster) continue;
213  else {
214  double cemc_ecore = cluster->get_ecore();
215  if(cemc_ecore<1.0) continue;
216  double cemc_x = cluster->get_x();
217  double cemc_y = cluster->get_y();
218  double cemc_z = cluster->get_z() - Zvtx;
219  double cemc_r = cluster->get_r();
220  double cemc_eta = asinh( cemc_z / cemc_r );
221  double cemc_phi = atan2( cemc_y, cemc_x );
222  double tmpdist = sqrt(pow((cemc_eta-track_eta),2)+pow((cemc_phi-track_phi),2));
223  if(tmpdist<dist) {
224  dist = tmpdist; returnCluster = cluster; dphi = fabs(cemc_phi-track_phi); deta = fabs(cemc_eta-track_eta);
225  }
226  }
227  }
228 
229  return returnCluster;
230 }
231 
232 //--------------------------------------------------------------------------------
233 
234 int FilterEventsUpsilon::process_event(PHCompositeNode *topNode) {
235 
236  EventNumber++;
237  int ngood = 0;
238  bool goodevent = false;
239 
240 cout << "getting nodes..." << endl;
241  GetNodes(topNode);
242 
243  GlobalVertexMap* _global_vtxmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
244  if(!_global_vtxmap) { cerr << PHWHERE << "ERROR: GlobalVertexMap node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
245 
246  double Zvtx = 0.;
247  for (GlobalVertexMap::Iter iter = _global_vtxmap->begin(); iter != _global_vtxmap->end(); ++iter)
248  {
249  GlobalVertex *vtx = iter->second;
250  if(vtx->get_id()==1) { Zvtx = vtx->get_z(); } // BBC vertex
251  }
252  cout << "Global BBC vertex Z = " << Zvtx << endl;
253 
254  std::vector<RawCluster *> goodclusters;
255  std::vector<TrkrCluster*> vclussilicon;
256  std::vector<TrkrCluster*> vclustpc;
257  std::vector<TrkrDefs::cluskey> vcluskeysilicon;
258  std::vector<TrkrDefs::cluskey> vcluskeytpc;
259  std::vector<SvtxTrackSeed_v1*> v_svtx_trackseed;
260  std::vector<TrackSeed*> v_silicon_trackseed;
261  std::vector<TrackSeed*> v_tpc_trackseed;
262 
263  cout << "Total number of tracks = " << _trackmap->size() << endl;
264  cout << "Total number of CEMC clusters = " << _cemc_clusters->size() << endl;
265  cout << "Total number of svtx track seeds: " << _trackseedcontainer_svtx->size() << endl;
266  cout << "Total number of silicon track seeds: " << _trackseedcontainer_silicon->size() << endl;
267  cout << "Total number of tpc track seeds: " << _trackseedcontainer_tpc->size() << endl;
268  cout << "Total number of TRKR clusters: " << _trkrclusters->size() << endl;
269 
270 // Start loop over tracks;
271  for (SvtxTrackMap::Iter iter = _trackmap->begin(); iter != _trackmap->end(); ++iter)
272  {
273  SvtxTrack *track = iter->second;
274 
275  double px = track->get_px();
276  double py = track->get_py();
277  double pz = track->get_pz();
278  double pt = sqrt(px * px + py * py);
279  if(pt<2.0) continue;
280  double mom = sqrt(px * px + py * py + pz * pz);
281 
282  double cemc_dphi = 99999.;
283  double cemc_deta = 99999.;
284  RawCluster* clus = MatchClusterCEMC(track, _cemc_clusters, cemc_dphi, cemc_deta, Zvtx);
285  if(!clus) continue;
286  double cemc_ecore = clus->get_ecore();
287  if(cemc_ecore/mom < 0.7) continue;
288 
289  ngood++;
290  goodevent = true;
291  SvtxTrack* tmp =_trackmap_ee->insert(track);
292  if(!tmp) cout << "ERROR: Failed to insert a track." << endl;
293 
294  cout << " Track: " << pt << endl;
295 
296  TrackSeed* trackseed_silicon = track->get_silicon_seed();
297  double trackseed_silicon_pt = trackseed_silicon->get_pt();
298  cout << " Silicon seed: " << trackseed_silicon << " " << trackseed_silicon_pt << " " << trackseed_silicon->size_cluster_keys() << endl;
299  TrackSeed_v1* tmpseedsilicon = (TrackSeed_v1*)trackseed_silicon->CloneMe();
300  v_silicon_trackseed.push_back(tmpseedsilicon);
301 
302  TrackSeed* trackseed_tpc = track->get_tpc_seed();
303  double trackseed_tpc_pt = trackseed_tpc->get_pt();
304  cout << " TPC seed: " << trackseed_tpc << " " << trackseed_tpc_pt << " " << trackseed_tpc->size_cluster_keys() << endl;
305  TrackSeed_v1* tmpseedtpc = (TrackSeed_v1*)trackseed_tpc->CloneMe();
306  v_tpc_trackseed.push_back(tmpseedtpc);
307 
308 // Find SVTX seed for this track
309  for(TrackSeedContainer::ConstIter seediter = _trackseedcontainer_svtx->begin(); seediter != _trackseedcontainer_svtx->end(); ++seediter)
310  {
311  SvtxTrackSeed_v1* seed = (SvtxTrackSeed_v1*)*seediter;
312  //TrackSeed* seed = *seediter;
313  bool foundsiliconseed = false;
314  bool foundtpcseed = false;
315  unsigned int siliconid = seed->get_silicon_seed_index();
316  unsigned int tpcid = seed->get_tpc_seed_index();
317  //cout << " SVTX seed: " << seed << " " << siliconid << " " << tpcid << endl;
318  TrackSeed* tmptpcseed = _trackseedcontainer_tpc->get(tpcid);
319  double tmptpcseed_pt = tmptpcseed->get_pt();
320  TrackSeed* tmpsiliconseed = _trackseedcontainer_silicon->get(siliconid);
321  double tmpsiliconseed_pt = tmpsiliconseed->get_pt();
322  if(tmpsiliconseed_pt == trackseed_silicon_pt) {foundsiliconseed = true; cout << " Found silicon seed " << tmpsiliconseed << endl;}
323  if(tmptpcseed_pt == trackseed_tpc_pt) {foundtpcseed = true; cout << " Found tpc seed " << tmptpcseed << endl;}
324  if(foundsiliconseed && foundtpcseed) {
325  cout << " THIS IS THE ONE: " << seed << " " << siliconid << " " << tpcid << endl;
326  SvtxTrackSeed_v1* tmpseed = (SvtxTrackSeed_v1*)seed->CloneMe();
327  v_svtx_trackseed.push_back(tmpseed);
328  }
329  }
330 
331 // Find all TRKR clusters for this track
332  for(auto clusiter = trackseed_tpc->begin_cluster_keys(); clusiter != trackseed_tpc->end_cluster_keys(); ++clusiter)
333  {
334  //auto key = *clusiter;
335  TrkrDefs::cluskey key = *clusiter;
336  TrkrCluster* clus = _trkrclusters->findCluster(key);
337  bool isinserted = false;
338  for(unsigned int i=0; i<vclustpc.size(); i++) {if(clus==vclustpc[i]) {isinserted=true; break;}}
339  if(!isinserted) {TrkrCluster* newclus = (TrkrClusterv4*)clus->CloneMe(); vclustpc.push_back(newclus); vcluskeytpc.push_back(key);}
340  }
341  for(auto clusiter = trackseed_silicon->begin_cluster_keys(); clusiter != trackseed_silicon->end_cluster_keys(); ++clusiter)
342  {
343  //auto key = *clusiter;
344  TrkrDefs::cluskey key = *clusiter;
345  TrkrCluster* clus = _trkrclusters->findCluster(key);
346  bool isinserted = false;
347  for(unsigned int i=0; i<vclussilicon.size(); i++) {if(clus==vclussilicon[i]) {isinserted=true; break;}}
348  if(!isinserted) {TrkrCluster* newclus = (TrkrClusterv4*)clus->CloneMe(); vclussilicon.push_back(newclus); vcluskeysilicon.push_back(key);}
349  }
350 
351  TVector3 proj = GetProjectionCEMC(track);
352  double track_x = proj(0);
353  double track_y = proj(1);
354  double track_z = proj(2) - Zvtx;
355  double track_r = sqrt(track_x*track_x+track_y*track_y);
356  double track_eta = asinh( track_z / track_r );
357  double track_phi = atan2( track_y, track_x );
358 
359 // Find all CEMC clusters around this track projection
360  RawClusterContainer::Range begin_end = _cemc_clusters->getClusters();
361  RawClusterContainer::Iterator clusiter;
362  for (clusiter = begin_end.first; clusiter != begin_end.second; ++clusiter)
363  {
364  RawCluster* cluster = clusiter->second;
365  if(!cluster) { cout << "ERROR: bad cluster pointer = " << cluster << endl; continue; }
366  else {
367  double cemc_ecore = cluster->get_ecore();
368  if(cemc_ecore<1.0) continue;
369  double cemc_x = cluster->get_x();
370  double cemc_y = cluster->get_y();
371  double cemc_z = cluster->get_z() - Zvtx;
372  double cemc_r = cluster->get_r();
373  double cemc_eta = asinh( cemc_z / cemc_r );
374  double cemc_phi = atan2( cemc_y, cemc_x );
375  double dist = sqrt(pow(cemc_phi-track_phi,2)+pow(cemc_eta-track_eta,2));
376  if(dist<0.1) {
377  bool isinserted = false;
378  for(unsigned int i=0; i<goodclusters.size(); i++) {if(cluster==goodclusters[i]) {isinserted=true; break;}}
379  if(!isinserted) {RawCluster* newcluster = (RawClusterv1*)cluster->CloneMe(); goodclusters.push_back(newcluster);}
380  }
381  }
382  } // end loop over cemc clusters
383 
384  } // end loop over tracks
385 
386  cout << " Number of CEMC clusters for output = " << goodclusters.size() << endl;
387  cout << " Number of svtx seeds for output = " << v_svtx_trackseed.size() << endl;
388  cout << " Number of tpc seeds for output = " << v_tpc_trackseed.size() << endl;
389  cout << " Number of silicon seeds for output = " << v_silicon_trackseed.size() << endl;
390 
391 // Fill in selected CEMC clusters
392  _cemc_clusters->Reset();
393  for (unsigned int cl = 0; cl < goodclusters.size(); cl++) { _cemc_clusters->AddCluster(goodclusters[cl]); }
394 
395 
396  cout << " Number of TRKR clusters for output: " << vclussilicon.size() << " " << vclustpc.size() << endl;
397  cout << " Number of TRKR cluster keys for output: " << vcluskeysilicon.size() << " " << vcluskeytpc.size() << endl;
398  _trkrclusters->Reset();
399  for(unsigned int cl = 0; cl < vclussilicon.size(); cl++)
400  { _trkrclusters->addClusterSpecifyKey(vcluskeysilicon[cl], vclussilicon[cl]); }
401  for(unsigned int cl = 0; cl < vclustpc.size(); cl++)
402  { _trkrclusters->addClusterSpecifyKey(vcluskeytpc[cl], vclustpc[cl]); }
403  cout << " New TRKR_CLUSTER size = " << _trkrclusters->size() << endl;
404 
405 
406  _trackseedcontainer_svtx->Reset();
407  for(unsigned int cl = 0; cl < v_svtx_trackseed.size(); cl++)
408  { _trackseedcontainer_svtx->insert(v_svtx_trackseed[cl]); }
409 
410 /*
411  _trackseedcontainer_silicon->Reset();
412  for(unsigned int cl = 0; cl < v_silicon_trackseed.size(); cl++)
413  { _trackseedcontainer_silicon->insert(v_silicon_trackseed[cl]); }
414  _trackseedcontainer_tpc->Reset();
415  for(unsigned int cl = 0; cl < v_tpc_trackseed.size(); cl++)
416  { _trackseedcontainer_tpc->insert(v_tpc_trackseed[cl]); }
417 */
418 
419  cout << " New numbers of seeds = " << _trackseedcontainer_svtx->size() << " " << _trackseedcontainer_silicon->size() << " " << _trackseedcontainer_tpc->size() << endl;
420  if(ngood>=2) { goodEventNumber++; }
421 
422  cout << ngood << " " << EventNumber << " " << goodEventNumber << endl;
423  if(goodevent) {return Fun4AllReturnCodes::EVENT_OK;} else {return Fun4AllReturnCodes::ABORTEVENT;}
424 }
425 
426 //---------------------------------------------------------------------------
427 
428 int FilterEventsUpsilon::End(PHCompositeNode *topNode)
429 {
430  cout << "Number of scanned events = " << EventNumber << endl;
431  cout << "Number of good events = " << goodEventNumber << endl;
432  return Fun4AllReturnCodes::EVENT_OK;
433 }
434