ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CreateCZHitContainer.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CreateCZHitContainer.cc
1 #include "CreateCZHitContainer.h"
2 
3 #include <g4main/PHG4Hit.h>
5 #include <g4main/PHG4Hitv1.h>
6 #include <g4main/PHG4Particle.h>
8 
11 #include <fun4all/SubsysReco.h> // for SubsysReco
12 
13 #include <TMath.h>
14 #include <TSystem.h>
15 #include <phool/PHCompositeNode.h>
16 #include <phool/PHIODataNode.h> // for PHIODataNode
17 #include <phool/PHNode.h> // for PHNode
18 #include <phool/PHNodeIterator.h> // for PHNodeIterator
19 #include <phool/PHObject.h> // for PHObject
20 #include <phool/getClass.h>
21 
22 #include <sstream>
23 #include <utility> // for pair
24 
25 using namespace std;
26 
28  : SubsysReco("CreateCZHitContainer" + name)
29  , _node_postfix(name)
30  , _truth_container(NULL)
31  , _hit_C(nullptr)
32  , _hit_Z(nullptr)
33  , _hit_CZ(nullptr)
34 {
35 }
36 
38 {
39 }
40 
42 {
43  _truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
44  if (!_truth_container)
45  {
46  cout << PHWHERE << " PHG4TruthInfoContainer node not found on node tree"
47  << endl;
49  }
50  ostringstream nodename;
51  nodename.str("");
52  nodename << "G4HIT_CZ" << _node_postfix;
53  PHNodeIterator iter(topNode);
54  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", _node_postfix));
55  PHG4HitContainer *cylinder_hits = findNode::getClass<PHG4HitContainer>(topNode, nodename.str());
56  if (!cylinder_hits)
57  {
58  dstNode->addNode(new PHIODataNode<PHObject>(cylinder_hits = new PHG4HitContainer(nodename.str()), nodename.str(), "PHObject"));
59  }
60  //cylinder_hits->AddLayer(GetLayer());
62 }
63 
65 {
66  ostringstream nodename;
67  nodename.str("");
68  nodename << "G4HIT_" << _node_postfix;
69  PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode, nodename.str().c_str());
70  nodename.str("");
71  nodename.clear();
72  nodename << "G4HIT_CZ" << _node_postfix;
73  PHG4HitContainer *cylinder_hits = findNode::getClass<PHG4HitContainer>(topNode, nodename.str().c_str());
74  if (!cylinder_hits)
75  {
76  cerr << "can't find " << nodename.str() << endl;
77  gSystem->Exit(1);
78  }
79  if (hits)
80  {
82  for (PHG4TruthInfoContainer::ConstIterator itr = itr_range.first;
83  itr != itr_range.second; ++itr)
84  {
85  PHG4Particle *particle = itr->second;
86  if (!particle) continue;
87  int trkid = particle->get_track_id();
88  for (PHG4HitContainer::LayerIter layerit =
89  hits->getLayers().first;
90  layerit != hits->getLayers().second; layerit++)
91  {
92  cylinder_hits->AddLayer(*layerit);
93  PHG4HitContainer::ConstRange hit_range = hits->getHits(*layerit);
94  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
95  {
96  if (hit_iter->second->get_trkid() != trkid) continue;
97  if (hit_iter->second->get_hit_type() == 1)
98  _hit_C = hit_iter->second;
99  else if (hit_iter->second->get_hit_type() == 2)
100  _hit_Z = hit_iter->second;
101  }
103  if (_hit_CZ)
104  {
105  cylinder_hits->AddHit(*layerit, _hit_CZ);
106  }
107  }
108  }
109  }
111 }
112 
114 {
115  if (h1 == nullptr || h2 == nullptr) return nullptr;
116  if (h1->get_trkid() != h2->get_trkid()) return nullptr;
117  if (h1->get_layer() != h2->get_layer()) return nullptr;
118  if (!(h1->get_hit_type() == 1 && h2->get_hit_type() == 2)) return nullptr;
119  if (!(h1->get_eion() > 0 && h2->get_eion() > 0)) return nullptr;
120  PHG4Hit *merged = new PHG4Hitv1();
121  merged->set_layer(h1->get_layer());
122 
123  float R_in = TMath::Sqrt(h1->get_x(0) * h1->get_x(0) + h1->get_y(0) * h1->get_y(0));
124  float R_out = TMath::Sqrt(h2->get_x(1) * h2->get_x(1) + h2->get_y(1) * h2->get_y(1));
125  float phi_measured = TMath::ATan2(h2->get_avg_y(), h2->get_avg_x());
126 
127  merged->set_x(0, R_in * TMath::Cos(phi_measured));
128  merged->set_y(0, R_in * TMath::Sin(phi_measured));
129  merged->set_z(0, h1->get_z(0));
130  merged->set_x(1, R_out * TMath::Cos(phi_measured));
131  merged->set_y(1, R_out * TMath::Sin(phi_measured));
132  merged->set_z(1, h1->get_z(1));
133  float t0 = h1->get_t(0);
134  if (h2->get_t(0) < t0) t0 = h2->get_t(0);
135  float t1 = h1->get_t(1);
136  if (h2->get_t(1) > t1) t1 = h2->get_t(1);
137  merged->set_t(0, t0);
138  merged->set_t(1, t1);
139  merged->set_trkid(h1->get_trkid());
140  return merged;
141 }