17 #include <phparameter/PHParameterContainerInterface.h>
51 cout <<
"Creating PHG4MvtxHitReco for name = " << name << endl;
63 std::cout <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
70 cout <<
"Could not locate g4 hit node " <<
hitnodename << endl;
78 cout <<
"Could not locate geometry node " <<
geonodename << endl;
86 auto hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
103 auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode,
"TRKR_HITTRUTHASSOC");
130 cout <<
"Could not locate g4 hit node " <<
hitnodename << endl;
137 cout <<
"Could not locate geometry node " <<
geonodename << endl;
142 auto trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
143 if (!trkrhitsetcontainer)
145 cout <<
"Could not locate TRKR_HITSET node, quit! " << endl;
150 auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode,
"TRKR_HITTRUTHASSOC");
153 cout <<
"Could not locate TRKR_HITTRUTHASSOC node, quit! " << endl;
159 pair<PHG4HitContainer::LayerIter, PHG4HitContainer::LayerIter> layer_begin_end = g4hit->
getLayers();
160 for (layer = layer_begin_end.first; layer != layer_begin_end.second; ++layer)
178 double xpixw_half = xpixw / 2.0;
180 double zpixw_half = zpixw / 2.0;
181 int maxNX = layergeom->
get_NX();
182 int maxNZ = layergeom->
get_NZ();
185 for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
189 hiter->second->print();
194 cout <<
PHWHERE <<
"Mvtx layers only go up to three! Quit." << endl;
198 cout <<
" layer " << *layer <<
" t0 " << hiter->second->get_t(0) <<
" t1 " << hiter->second->get_t(1)
201 if (hiter->second->get_t(0) >
tmin_max[*
layer].second)
continue;
202 if (hiter->second->get_t(1) <
tmin_max[*
layer].first)
continue;
210 TVector3 local_in(hiter->second->get_local_x(0), hiter->second->get_local_y(0), hiter->second->get_local_z(0));
211 TVector3 local_out(hiter->second->get_local_x(1), hiter->second->get_local_y(1), hiter->second->get_local_z(1));
212 TVector3 midpoint((local_in.X() + local_out.X()) / 2.0, (local_in.Y() + local_out.Y()) / 2.0, (local_in.Z() + local_out.Z()) / 2.0);
217 <<
" world entry point position: " << hiter->second->get_x(0) <<
" " << hiter->second->get_y(0) <<
" " << hiter->second->get_z(0) << endl;
218 cout <<
" world exit point position: " << hiter->second->get_x(1) <<
" " << hiter->second->get_y(1) <<
" " << hiter->second->get_z(1) << endl;
219 cout <<
" local coords of entry point from G4 " << hiter->second->get_local_x(0) <<
" " << hiter->second->get_local_y(0) <<
" " << hiter->second->get_local_z(0) << endl;
220 TVector3 world_in(hiter->second->get_x(0), hiter->second->get_y(0), hiter->second->get_z(0));
221 TVector3 local_in_check = layergeom->
get_local_from_world_coords(stave_number, half_stave_number, module_number, chip_number, world_in);
222 cout <<
" local coords of entry point from geom (a check) " << local_in_check.X() <<
" " << local_in_check.Y() <<
" " << local_in_check.Z() << endl;
223 cout <<
" local coords of exit point from G4 " << hiter->second->get_local_x(1) <<
" " << hiter->second->get_local_y(1) <<
" " << hiter->second->get_local_z(1) << endl;
224 cout <<
" local coords of exit point from geom (a check) " << local_out.X() <<
" " << local_out.Y() <<
" " << local_out.Z() << endl;
231 TVector3 location_in = layergeom->
get_world_from_local_coords(stave_number, half_stave_number, module_number, chip_number, local_in);
232 TVector3 location_out = layergeom->
get_world_from_local_coords(stave_number, half_stave_number, module_number, chip_number, local_out);
235 <<
" PHG4MvtxHitReco: Found world entry location from geometry for "
236 <<
" stave number " << stave_number
237 <<
" half stave number " << half_stave_number
238 <<
" module number " << module_number
239 <<
" chip number " << chip_number
241 <<
" x = " << location_in.X()
242 <<
" y = " << location_in.Y()
243 <<
" z = " << location_in.Z()
244 <<
" radius " << sqrt(pow(location_in.X(), 2) + pow(location_in.Y(), 2))
245 <<
" angle " << atan(location_in.Y() / location_in.X())
247 cout <<
" PHG4MvtxHitReco: The world entry location from G4 was "
249 <<
" x = " << hiter->second->get_x(0)
250 <<
" y = " << hiter->second->get_y(0)
251 <<
" z = " << hiter->second->get_z(0)
252 <<
" radius " << sqrt(pow(hiter->second->get_x(0), 2) + pow(hiter->second->get_y(0), 2))
253 <<
" angle " << atan(hiter->second->get_y(0) / hiter->second->get_x(0))
255 cout <<
" difference in x = " << hiter->second->get_x(0) - location_in.X()
256 <<
" difference in y = " << hiter->second->get_y(0) - location_in.Y()
257 <<
" difference in z = " << hiter->second->get_z(0) - location_in.Z()
258 <<
" difference in radius = " << sqrt(pow(hiter->second->get_x(0), 2) + pow(hiter->second->get_y(0), 2)) - sqrt(pow(location_in.X(), 2) + pow(location_in.Y(), 2))
259 <<
" in angle = " << atan(hiter->second->get_y(0) / hiter->second->get_x(0)) - atan(location_in.Y() / location_in.X())
263 cout <<
" PHG4MvtxHitReco: Found world exit location from geometry for "
264 <<
" stave number " << stave_number
265 <<
" half stave number " << half_stave_number
266 <<
" module number" << module_number
268 <<
" x = " << location_out.X()
269 <<
" y = " << location_out.Y()
270 <<
" z = " << location_out.Z()
271 <<
" radius " << sqrt(pow(location_out.X(), 2) + pow(location_out.Y(), 2))
272 <<
" angle " << atan(location_out.Y() / location_out.X())
274 cout <<
" PHG4MvtxHitReco: The world exit location from G4 was "
276 <<
" x = " << hiter->second->get_x(1)
277 <<
" y = " << hiter->second->get_y(1)
278 <<
" z = " << hiter->second->get_z(1)
279 <<
" radius " << sqrt(pow(hiter->second->get_x(1), 2) + pow(hiter->second->get_y(1), 2))
280 <<
" angle " << atan(hiter->second->get_y(1) / hiter->second->get_x(1))
282 cout <<
" difference in radius = " << sqrt(pow(hiter->second->get_x(1), 2) + pow(hiter->second->get_y(1), 2)) - sqrt(pow(location_out.X(), 2) + pow(location_out.Y(), 2))
283 <<
" in angle = " << atan(hiter->second->get_y(1) / hiter->second->get_x(1)) - atan(location_out.Y() / location_out.X())
293 if (pixel_number_in < 0 || pixel_number_out < 0)
295 cout <<
"Oops! got negative pixel number in layer " << layergeom->
get_layer()
296 <<
" pixel_number_in " << pixel_number_in
297 <<
" pixel_number_out " << pixel_number_out
298 <<
" local_in = " << local_in.X() <<
" " << local_in.Y() <<
" " << local_in.Z()
299 <<
" local_out = " << local_out.X() <<
" " << local_out.Y() <<
" " << local_out.Z()
304 cout <<
"entry pixel number " << pixel_number_in <<
" exit pixel number " << pixel_number_out << endl;
309 vector<pair<double, double> > venergy;
327 TVector3 pathvec = local_in - local_out;
335 double diffusion_width_max = 25.0e-04;
336 double diffusion_width_min = 8.0e-04;
338 double ydrift_max = pathvec.Y();
350 int xbin_max, xbin_min;
352 if (xbin_in > xbin_out)
354 xbin_max = xbin_in + nadd;
355 xbin_min = xbin_out - nadd;
359 xbin_max = xbin_out + nadd;
360 xbin_min = xbin_in - nadd;
363 int zbin_max, zbin_min;
364 if (zbin_in > zbin_out)
366 zbin_max = zbin_in + nadd;
367 zbin_min = zbin_out - nadd;
371 zbin_max = zbin_out + nadd;
372 zbin_min = zbin_in - nadd;
377 if (xbin_min < 0) xbin_min = 0;
378 if (zbin_min < 0) zbin_min = 0;
379 if (xbin_max >= maxNX) xbin_max = maxNX-1;
380 if (zbin_max >= maxNZ) xbin_max = maxNZ-1;
384 cout <<
" xbin_in " << xbin_in <<
" xbin_out " << xbin_out <<
" xbin_min " << xbin_min <<
" xbin_max " << xbin_max << endl;
385 cout <<
" zbin_in " << zbin_in <<
" zbin_out " << zbin_out <<
" zbin_min " << zbin_min <<
" zbin_max " << zbin_max << endl;
390 if (xbin_max - xbin_min > 12 || zbin_max - zbin_min > 12)
394 double pixenergy[12][12] = {};
395 double pixeion[12][12] = {};
398 for (
int i = 0; i < nsegments; i++)
403 double interval = 2 * (double) i + 1;
404 double frac = interval / (double) (2 * nsegments);
405 TVector3 segvec(pathvec.X() * frac, pathvec.Y() * frac, pathvec.Z() * frac);
406 segvec = segvec + local_out;
410 double ydrift = segvec.Y() - local_out.Y();
414 double ydiffusion_radius = diffusion_width_min + (ydrift / ydrift_max) * (diffusion_width_max - diffusion_width_min);
417 cout <<
" segment " << i
418 <<
" interval " << interval
420 <<
" local_in.X " << local_in.X()
421 <<
" local_in.Z " << local_in.Z()
422 <<
" local_in.Y " << local_in.Y()
423 <<
" pathvec.X " << pathvec.X()
424 <<
" pathvec.Z " << pathvec.Z()
425 <<
" pathvec.Y " << pathvec.Y()
426 <<
" segvec.X " << segvec.X()
427 <<
" segvec.Z " << segvec.Z()
428 <<
" segvec.Y " << segvec.Y()
429 <<
" ydrift " << ydrift
430 <<
" ydrift_max " << ydrift_max
431 <<
" ydiffusion_radius " << ydiffusion_radius
435 for (
int ix = xbin_min; ix <= xbin_max; ix++)
437 for (
int iz = zbin_min; iz <= zbin_max; iz++)
444 cout <<
" pixnum < 0 , pixnum = " << pixnum << endl;
445 cout <<
" ix " << ix <<
" iz " << iz << endl;
446 cout <<
" xbin_min " << xbin_min <<
" zbin_min " << zbin_min
447 <<
" xbin_max " << xbin_max <<
" zbin_max " << zbin_max
449 cout <<
" maxNX " << maxNX <<
" maxNZ " << maxNZ
455 double x1 = tmp.X() - xpixw_half;
456 double z1 = tmp.Z() + zpixw_half;
457 double x2 = tmp.X() + xpixw_half;
458 double z2 = tmp.Z() - zpixw_half;
464 pixenergy[ix - xbin_min][iz - zbin_min] += pixarea_frac * hiter->second->get_edep() / (float) nsegments;
467 pixeion[ix - xbin_min][iz - zbin_min] += pixarea_frac * hiter->second->get_eion() / (float) nsegments;
471 cout <<
" pixnum " << pixnum <<
" xbin " << ix <<
" zbin " << iz
472 <<
" pixel_area fraction of circle " << pixarea_frac <<
" accumulated pixel energy " << pixenergy[ix - xbin_min][iz - zbin_min]
480 for (
int ix = xbin_min; ix <= xbin_max; ix++)
482 for (
int iz = zbin_min; iz <= zbin_max; iz++)
484 if (pixenergy[ix - xbin_min][iz - zbin_min] > 0.0)
487 vpixel.push_back(pixnum);
490 pair<double, double> tmppair = make_pair(pixenergy[ix - xbin_min][iz - zbin_min], pixeion[ix - xbin_min][iz - zbin_min]);
491 venergy.push_back(tmppair);
493 cout <<
" Added pixel number " << pixnum <<
" xbin " << ix <<
" zbin " << iz <<
" to vectors with energy " << pixenergy[ix - xbin_min][iz - zbin_min] << endl;
503 for (
unsigned int i1 = 0; i1 < vpixel.size(); i1++)
517 hit = hitsetit->second->getHit(hitkey);
522 hitsetit->second->addHitSpecificKey(hitkey, hit);
526 hit->
addEnergy(venergy[i1].first * TrkrDefs::MvtxEnergyScaleup);
534 hittruthassoc->addAssoc(hitsetkey, hitkey, hiter->first);
544 cout <<
"From PHG4MvtxHitReco: " << endl;
545 hittruthassoc->identify();
554 std::map<int, std::pair<double, double> >::iterator
it =
tmin_max.find(detid);
557 tmin_max.insert(std::make_pair(detid, std::make_pair(tmin, tmax)));
559 cout <<
"PHG4MvtxHitReco: Set Mvtx timing window parameters from macro for layer = " << detid <<
" to tmin = " <<
tmin_max[detid].first <<
" tmax = " <<
tmin_max[detid].second << endl;
566 cout <<
"PHG4MvtxHitReco: Setting Mvtx timing window defaults to tmin = -5000 and tmax = 5000 ns" << endl;
567 for (
int ilayer = 0; ilayer < 3; ilayer++)
569 tmin_max.insert(std::make_pair(ilayer, std::make_pair(-5000, 5000)));