ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4TTLSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4TTLSteppingAction.cc
2 #include "PHG4TTLDetector.h"
3 
4 
6 class PHCompositeNode;
7 //____________________________________________________________________________..
9  : PHG4SteppingAction(detector->GetName())
10  , detector_(detector)
11  , hits_(nullptr)
12  , hit(nullptr)
13  , saveshower(nullptr)
14  , layer_id(-1)
15  , _isFwd_TTL(false)
16  , _N_phi_modules(0)
17  , _z_pos_TTL(0)
18  , _sensor_resolution_x(0)
19  , _sensor_resolution_y(0)
20 {
21 }
22 
24 {
25  // if the last hit was a zero energie deposit hit, it is just reset
26  // and the memory is still allocated, so we need to delete it here
27  // if the last hit was saved, hit is a nullptr pointer which are
28  // legal to delete (it results in a no operation)
29  delete hit;
30 }
31 
32 //____________________________________________________________________________..
34 {
35  // get volume of the current step
38 
40  // collect energy and track length step by step
42  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
43 
44  const G4Track* aTrack = aStep->GetTrack();
45 
46  TVector3 sensorPosition;
47  // make sure we are in a volume
48  if (detector_->IsInSectorActive(volume))
49  {
50  bool geantino = false;
51 
52  int moduleID = -1;
53  int layer = -1;
54  int sensor0 = -1;
55  int sensor1 = -1;
56  int idx_j = -1;
57  int idx_k = -1;
58  // the check for the pdg code speeds things up, I do not want to make
59  // an expensive string compare for every track when we know
60  // geantino or chargedgeantino has pid=0
61  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 && aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
62  {
63  geantino = true;
64  }
65  G4StepPoint* prePoint = aStep->GetPreStepPoint();
66  G4StepPoint* postPoint = aStep->GetPostStepPoint();
67  // cout << "track id " << aTrack->GetTrackID() << endl;
68  // cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
69  // cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
70  //layer_id is sector number
71  switch (prePoint->GetStepStatus())
72  {
73  case fGeomBoundary:
74  case fUndefined:
75  if (!hit)
76  {
77  hit = new PHG4Hitv1();
78  }
79  //here we set the entrance values in cm
80  hit->set_x(0, prePoint->GetPosition().x() / cm);
81  hit->set_y(0, prePoint->GetPosition().y() / cm);
82  hit->set_z(0, prePoint->GetPosition().z() / cm);
83  // time in ns
84  hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
85  //set the track ID
86  hit->set_trkid(aTrack->GetTrackID());
88  {
89  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
90  {
91  hit->set_trkid(pp->GetUserTrackId());
92  hit->set_shower_id(pp->GetShower()->get_id());
93  saveshower = pp->GetShower();
94  }
95  }
96 
97  // std::cout << std::endl;
98  CalculateSensorHitIndices(prePoint, moduleID,layer,sensor0,sensor1,idx_j, idx_k,sensorPosition);
99  hit->set_index_i(moduleID);
100  hit->set_index_j(layer);
101  hit->set_index_k(sensor0);
102  hit->set_index_l(sensor1);
103  hit->set_strip_z_index(idx_j);
104  hit->set_strip_y_index(idx_k);
105 
106  hit->set_local_x(0, sensorPosition.X());
107  hit->set_local_y(0, sensorPosition.Y());
108  hit->set_local_z(0, sensorPosition.Z());
109 
110  //set the initial energy deposit
111  hit->set_edep(0);
112  hit->set_eion(0); // only implemented for v5 otherwise empty
113  layer_id = 0;//aStep->GetPreStepPoint()->GetTouchable()->GetReplicaNumber(2);
114  // std::cout << "layerid: " << layer_id << std::endl;
115  // hit->set_light_yield(0);
116 
117  break;
118  default:
119  break;
120  }
121  // here we just update the exit values, it will be overwritten
122  // for every step until we leave the volume or the particle
123  // ceases to exist
124  hit->set_x(1, postPoint->GetPosition().x() / cm);
125  hit->set_y(1, postPoint->GetPosition().y() / cm);
126  hit->set_z(1, postPoint->GetPosition().z() / cm);
127 
128  hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
129  //sum up the energy to get total deposited
130  hit->set_edep(hit->get_edep() + edep);
131  // std::cout << "energy: " << hit->get_edep() + edep << std::endl;
132  hit->set_eion(hit->get_eion() + eion);
133  hit->set_path_length(aTrack->GetTrackLength() / cm);
134  if (geantino)
135  {
136  hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
137  }
138  if (edep > 0)
139  {
141  {
142  if (PHG4TrackUserInfoV1* pp =
143  dynamic_cast<PHG4TrackUserInfoV1*>(p))
144  {
145  pp->SetKeep(1); // we want to keep the track
146  }
147  }
148  }
149  // if any of these conditions is true this is the last step in
150  // this volume and we need to save the hit
151  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
152  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
153  // (not sure if this will ever be the case)
154  // aTrack->GetTrackStatus() == fStopAndKill: track ends
155  if (postPoint->GetStepStatus() == fGeomBoundary ||
156  postPoint->GetStepStatus() == fWorldBoundary ||
157  postPoint->GetStepStatus() == fAtRestDoItProc ||
158  aTrack->GetTrackStatus() == fStopAndKill)
159  {
160  // save only hits with energy deposit (or -1 for geantino)
161  if (hit->get_edep())
162  {
164  if (saveshower)
165  {
167  }
168  // ownership has been transferred to container, set to null
169  // so we will create a new hit for the next track
170  hit = nullptr;
171  }
172  else
173  {
174  // if this hit has no energy deposit, just reset it for reuse
175  // this means we have to delete it in the dtor. If this was
176  // the last hit we processed the memory is still allocated
177  hit->Reset();
178  }
179  }
180 
181  // hit->identify();
182  // return true to indicate the hit was used
183  return true;
184  }
185  else
186  {
187  return false;
188  }
189 }
190 
191 
192 void PHG4TTLSteppingAction::CalculateSensorHitIndices(G4StepPoint* prePoint, int& module_ret, int& layer, int& sensor_0, int& sensor_1, int& j, int& k, TVector3 &sensorposition)
193 {
194  int module_ID = -1; //The j and k indices for the scintillator / tower
195  int layer_ID = -1; //The j and k indices for the scintillator / tower
196  int sensor_ID_0 = -1; //The j and k indices for the scintillator / tower
197  int sensor_ID_1 = -1; //The j and k indices for the scintillator / tower
198  int hit_j_0 = 0; //The j and k indices for the scintillator / tower
199  int hit_k_0 = 0; //The j and k indices for the scintillator / tower
200 
201 
202  G4double baseplate_length = 43.1 * mm;
203  G4double baseplate_width = 56.5 * mm / 2;
204  G4double baseSH_width = baseplate_width / 2;
205  G4double _module_x_dimension = baseplate_length;
206  G4double _module_y_dimension = baseplate_width + baseSH_width;
207 
208  G4double sensor_y_dimension = 21.2 * mm;
209  G4double sensor_x_dimension = 42.0 * mm;
210 
211  _sensor_resolution_x = (500e-4 / sqrt(12)) * cm; // in mm
212  _sensor_resolution_y = (500e-4 / sqrt(12)) * cm; // in mm
213 
214  TVector3 prePointVec(prePoint->GetPosition().x(),prePoint->GetPosition().y(),prePoint->GetPosition().z());
215  if(_N_phi_modules>0){
216  float prePoint_Phi = prePointVec.Phi()+M_PI;
217  module_ID = (int) prePoint_Phi/(2*M_PI/_N_phi_modules);
218  }
219  if(prePoint->GetPosition().z()>0){
220  layer_ID = prePoint->GetPosition().z()>_z_pos_TTL ? 1 : 0;
221  } else {
222  layer_ID = prePoint->GetPosition().z()<_z_pos_TTL ? 1 : 0;
223  }
224  if(_isFwd_TTL){
225  // front face starts with sensors
226  sensor_ID_0 = (int) ( ( ( prePoint->GetPosition().x() + (prePoint->GetPosition().x()<0 ? -_module_x_dimension /2 : _module_x_dimension/2) ) ) / _module_x_dimension );
227  sensor_ID_1 = (int) ( ( ( prePoint->GetPosition().y() + (prePoint->GetPosition().y()<0 ? -_module_y_dimension /2 : _module_y_dimension/2) ) ) / _module_y_dimension );
228  // calculate x and y position of bottom left corner of sensor
229  float sensorcorner_x = sensor_ID_0*_module_x_dimension - sensor_x_dimension/2;
230  float sensorcorner_y = 0;
231  if((sensor_ID_1%2==0 && layer_ID==0) || (sensor_ID_1%2!=0 && layer_ID==1) ){
232  // even sensor counts are located at the bottom of the module
233  sensorcorner_y = sensor_ID_1*_module_y_dimension - _module_y_dimension/2 + (0.1 * mm / 2);
234  } else {
235  // odd sensor counts are located at the top of the module
236  sensorcorner_y = sensor_ID_1*_module_y_dimension + _module_y_dimension/2 - sensor_y_dimension - (0.1 * mm / 2);
237  }
238  // std::cout << "\tcorner_x " << sensorcorner_x << "\tposition_x " << prePoint->GetPosition().x() << "\treso " << _sensor_resolution_x << std::endl;
239  // std::cout << "\tcorner_y " << sensorcorner_y << "\tposition_y " << prePoint->GetPosition().y() << "\treso " << _sensor_resolution_y << std::endl;
240  hit_j_0 = (int) ( ( prePoint->GetPosition().x() - sensorcorner_x) / _sensor_resolution_x );
241  sensorposition.SetX( (hit_j_0 * _sensor_resolution_x) + sensorcorner_x );
242  hit_k_0 = (int) ( ( prePoint->GetPosition().y() - sensorcorner_y) / _sensor_resolution_y );
243  sensorposition.SetY( (hit_k_0 * _sensor_resolution_y) + sensorcorner_y );
244 
245  sensorposition.SetZ( prePoint->GetPosition().z() );
246  }
247 
248 
249  module_ret = module_ID;
250  layer = layer_ID;
251  sensor_0 = sensor_ID_0;
252  sensor_1 = sensor_ID_1;
253  j = hit_j_0;
254  k = hit_k_0;
255  // if(_isFwd_TTL){
256  // std::cout << "module " << module_ID << "\tlayer " << layer_ID << "\tsensor0 " << sensor_ID_0 << "\tsensor1 " << sensor_ID_1 << "\tj " << j << "\tk " << k << std::endl;
257  // }
258  return;
259 }
260 
261 
262 //____________________________________________________________________________..
264 {
265  std::string hitnodename;
266  if (detector_->SuperDetector() != "NONE")
267  {
268  hitnodename = "G4HIT_" + detector_->SuperDetector();
269  }
270  else
271  {
272  hitnodename = "G4HIT_" + detector_->GetName();
273  }
274 
275  //now look for the map and grab a pointer to it.
276  hits_ = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
277 
278  // if we do not find the node we need to make it.
279  if (!hits_)
280  {
281  std::cout << "PHG4TTLSteppingAction::SetTopNode - unable to find "
282  << hitnodename << std::endl;
283  }
284 }