17 #include <phparameter/PHParameterInterface.h>
81 std::cout <<
PHWHERE <<
"DST Node missing, exiting." << std::endl;
90 cout <<
Name() <<
"RUN Node missing, exiting." << endl;
97 cout <<
Name() <<
"PAR Node missing, exiting." << endl;
101 string paramnodename =
"G4CELLPARAM_" +
m_Detector;
116 std::cout <<
"Could not locate g4 hit node " <<
m_HitNodeName << std::endl;
120 auto hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
121 if (!hitsetcontainer)
137 auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode,
"TRKR_HITTRUTHASSOC");
157 std::cout <<
"Could not locate geometry node " <<
m_GeoNodeName << std::endl;
188 std::cout <<
"Could not locate g4 hit node " <<
m_HitNodeName << std::endl;
193 auto hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
194 if (!hitsetcontainer)
196 cout <<
"Could not locate TRKR_HITSET node, quit! " << endl;
201 auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode,
"TRKR_HITTRUTHASSOC");
204 cout <<
"Could not locate TRKR_HITTRUTHASSOC node, quit! " << endl;
211 std::cout <<
"Could not locate geometry node " <<
m_GeoNodeName << std::endl;
216 if (
Verbosity() > 2) cout <<
" PHG4InttHitReco: Loop over hits" << endl;
220 const int sphxlayer = hiter->second->get_detid();
226 if (hiter->second->get_t(0) >
m_Tmax)
228 if (hiter->second->get_t(1) <
m_Tmin)
232 double diffusion_width = 5.0e-04;
234 const int ladder_z_index = hiter->second->get_ladder_z_index();
235 const int ladder_phi_index = hiter->second->get_ladder_phi_index();
240 int strip_y_index_in, strip_z_index_in, strip_y_index_out, strip_z_index_out;
241 layergeom->
find_strip_index_values(ladder_z_index, hiter->second->get_local_y(0), hiter->second->get_local_z(0), strip_y_index_in, strip_z_index_in);
242 layergeom->
find_strip_index_values(ladder_z_index, hiter->second->get_local_y(1), hiter->second->get_local_z(1), strip_y_index_out, strip_z_index_out);
247 double check_location[3] = {-1, -1, -1};
249 cout <<
" G4 entry location = " << hiter->second->get_local_x(0) <<
" " << hiter->second->get_local_y(0) <<
" " << hiter->second->get_local_z(0) << endl;
250 cout <<
" Check entry location = " << check_location[0] <<
" " << check_location[1] <<
" " << check_location[2] << endl;
252 cout <<
" G4 exit location = " << hiter->second->get_local_x(1) <<
" " << hiter->second->get_local_y(1) <<
" " << hiter->second->get_local_z(1) << endl;
253 cout <<
" Check exit location = " << check_location[0] <<
" " << check_location[1] <<
" " << check_location[2] << endl;
257 int minstrip_z = strip_z_index_in;
258 int maxstrip_z = strip_z_index_out;
259 if (minstrip_z > maxstrip_z)
swap(minstrip_z, maxstrip_z);
261 int minstrip_y = strip_y_index_in;
262 int maxstrip_y = strip_y_index_out;
263 if (minstrip_y > maxstrip_y)
swap(minstrip_y, maxstrip_y);
270 vector<pair<double, double> > venergy;
284 if (maxstrip_y - minstrip_y > 12 || maxstrip_z - minstrip_z > 12)
289 double stripenergy[13][13] = {};
290 double stripeion[13][13] = {};
295 gsl_vector_set(
m_PathVec, 0, hiter->second->get_local_x(0));
296 gsl_vector_set(
m_PathVec, 1, hiter->second->get_local_y(0));
297 gsl_vector_set(
m_PathVec, 2, hiter->second->get_local_z(0));
298 gsl_vector_set(
m_LocalOutVec, 0, hiter->second->get_local_x(1));
299 gsl_vector_set(
m_LocalOutVec, 1, hiter->second->get_local_y(1));
300 gsl_vector_set(
m_LocalOutVec, 2, hiter->second->get_local_z(1));
302 for (
int i = 0; i < nsegments; i++)
307 double interval = 2 * (double) i + 1;
308 double frac = interval / (double) (2 * nsegments);
314 double diffusion_radius = diffusion_width;
317 cout <<
" segment " << i
318 <<
" interval " << interval
320 <<
" local_in.X " << hiter->second->get_local_x(0)
321 <<
" local_in.Z " << hiter->second->get_local_z(0)
322 <<
" local_in.Y " << hiter->second->get_local_y(0)
323 <<
" pathvec.X " << gsl_vector_get(
m_PathVec, 0)
324 <<
" pathvec.Z " << gsl_vector_get(
m_PathVec, 2)
325 <<
" pathvec.Y " << gsl_vector_get(
m_PathVec, 1)
328 <<
" segvec.Y " << gsl_vector_get(
m_SegmentVec, 1) << endl
329 <<
" diffusion_radius " << diffusion_radius
333 for (
int iz = minstrip_z; iz <= maxstrip_z; iz++)
335 for (
int iy = minstrip_y; iy <= maxstrip_y; iy++)
338 double location[3] = {-1, -1, -1};
350 stripenergy[iy - minstrip_y][iz - minstrip_z] += striparea_frac * hiter->second->get_edep() / (float) nsegments;
353 stripeion[iy - minstrip_y][iz - minstrip_z] += striparea_frac * hiter->second->get_eion() / (float) nsegments;
357 cout <<
" strip y index " << iy <<
" strip z index " << iz
358 <<
" strip area fraction of circle " << striparea_frac <<
" accumulated pixel energy " << stripenergy[iy - minstrip_y][iz - minstrip_z]
366 for (
int iz = minstrip_z; iz <= maxstrip_z; iz++)
368 for (
int iy = minstrip_y; iy <= maxstrip_y; iy++)
370 if (stripenergy[iy - minstrip_y][iz - minstrip_z] > 0.0)
374 pair<double, double> tmppair = make_pair(stripenergy[iy - minstrip_y][iz - minstrip_z], stripeion[iy - minstrip_y][iz - minstrip_z]);
375 venergy.push_back(tmppair);
377 cout <<
" Added ybin " << iy <<
" zbin " << iz <<
" to vectors with energy " << stripenergy[iy - minstrip_y][iz - minstrip_z] << endl;
386 for (
unsigned int i1 = 0; i1 < vybin.size(); i1++)
399 TrkrHit *hit = hitsetit->second->getHit(hitkey);
405 hitsetit->second->addHitSpecificKey(hitkey, hit);
410 cout <<
"add energy " << venergy[i1].first <<
" to intthit " << endl;
411 hit->
addEnergy(venergy[i1].first * TrkrDefs::InttEnergyScaleup);
414 hittruthassoc->addAssoc(hitsetkey, hitkey, hiter->first);
417 cout <<
"PHG4InttHitReco: added hit wirh hitsetkey " << hitsetkey <<
" hitkey " << hitkey <<
" g4hitkey " << hiter->first <<
" energy " << hit->
getEnergy() << endl;
425 cout <<
"From PHG4InttHitReco: " << endl;
426 hitsetcontainer->identify();
427 hittruthassoc->identify();