24 #include <phparameter/PHParameterInterface.h>
41 #include <gsl/gsl_randist.h>
42 #include <gsl/gsl_rng.h>
55 inline constexpr
T square(
const T&
x ) {
return x*
x; }
59 inline T gaus(
const T&
x,
const T& sigma )
60 {
return std::exp( -
square(x/sigma)/2 )/(sigma*std::sqrt(2*
M_PI)); }
64 inline T bind_angle(
const T&
angle )
66 if( angle >=
M_PI )
return angle - 2*
M_PI;
67 else if( angle < -M_PI )
return angle + 2*
M_PI;
73 inline T get_rectangular_fraction(
const T& xloc,
const T& sigma,
const T& pitch )
74 {
return (std::erf( (xloc + pitch/2)/(M_SQRT2*sigma) ) - std::erf( (xloc - pitch/2)/(M_SQRT2*sigma) ))/2; }
81 inline T get_zigzag_fraction(
const T& xloc,
const T& sigma,
const T& pitch )
85 (pitch - xloc)*(std::erf(xloc/(M_SQRT2*sigma)) - std::erf((xloc-pitch)/(M_SQRT2*sigma)))/(pitch*2)
86 + (gaus(xloc-pitch, sigma) - gaus(xloc, sigma))*
square(sigma)/pitch
89 + (pitch + xloc)*(std::erf((xloc+pitch)/(M_SQRT2*sigma)) - std::erf(xloc/(M_SQRT2*sigma)))/(pitch*2)
90 + (gaus(xloc+pitch, sigma) - gaus(xloc, sigma))*
square(sigma)/pitch;
99 , m_detector(detector)
103 m_rng.reset( gsl_rng_alloc(gsl_rng_mt19937) );
128 <<
"PHG4MicromegasHitReco::InitRun\n"
129 <<
" m_tmin: " <<
m_tmin <<
"ns, m_tmax: " <<
m_tmax <<
"ns\n"
131 <<
" m_gain: " <<
m_gain <<
"\n"
146 std::cout <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
151 auto hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
152 if (!hitsetcontainer)
154 std::cout <<
PHWHERE <<
"creating TRKR_HITSET." << std::endl;
162 dstNode->addNode(trkrnode);
168 trkrnode->addNode(newNode);
172 auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode,
"TRKR_HITTRUTHASSOC");
175 std::cout <<
PHWHERE <<
"creating TRKR_HITTRUTHASSOC." << std::endl;
183 dstNode->addNode(trkrnode);
188 trkrnode->addNode(newNode);
199 const std::string g4hitnodename =
"G4HIT_" +
m_detector;
200 auto g4hitcontainer = findNode::getClass<PHG4HitContainer>(topNode, g4hitnodename);
201 assert(g4hitcontainer);
205 auto geonode = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename.c_str());
209 auto trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
210 assert(trkrhitsetcontainer);
213 auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode,
"TRKR_HITTRUTHASSOC");
214 assert(hittruthassoc);
217 auto layer_range = g4hitcontainer->getLayers();
218 for(
auto layer_it = layer_range.first; layer_it != layer_range.second; ++layer_it )
221 const auto layer = *layer_it;
234 const auto mesh_local_y = layergeom->get_drift_direction() == MicromegasDefs::DriftDirection::OUTWARD ?
235 layergeom->get_thickness()/2:
236 -layergeom->get_thickness()/2;
251 for(
auto g4hit_it = g4hit_range.first; g4hit_it != g4hit_range.second; ++g4hit_it )
255 PHG4Hit* g4hit = g4hit_it->second;
262 TVector3 world_in( g4hit->
get_x(0), g4hit->
get_y(0), g4hit->
get_z(0) );
263 TVector3 world_out( g4hit->
get_x(1), g4hit->
get_y(1), g4hit->
get_z(1) );
270 const int tileid = layergeom->find_tile_cylindrical( (world_in+world_out)*0.5 );
271 if( tileid < 0 )
continue;
288 layergeom->convert_to_planar( tileid, &g4hit_copy );
289 world_in.SetXYZ( g4hit_copy.
get_x(0), g4hit_copy.
get_y(0), g4hit_copy.
get_z(0) );
290 world_out.SetXYZ( g4hit_copy.
get_x(1), g4hit_copy.
get_y(1), g4hit_copy.
get_z(1) );
292 const auto local_in = layergeom->get_local_from_world_coords( tileid, world_in );
293 const auto local_out = layergeom->get_local_from_world_coords( tileid, world_out );
297 if( !nprimary )
continue;
301 const auto hitset_it = trkrhitsetcontainer->findOrAddHitSet(hitsetkey);
304 using charge_map_t = std::map<int,double>;
305 charge_map_t total_charges;
308 for( uint ie = 0; ie < nprimary; ++ie )
312 const auto t = gsl_ran_flat(
m_rng.get(), 0.0, 1.0);
313 auto local = local_in*
t + local_out*(1.0-
t);
319 const double y =
local.y();
320 const double drift_distance =
std::abs(y - mesh_local_y);
322 const double diffusion_angle = gsl_ran_flat(
m_rng.get(), -
M_PI,
M_PI);
325 local += TVector3( diffusion*std::cos(diffusion_angle), 0, diffusion*std::sin(diffusion_angle) );
334 const auto sum = std::accumulate( fractions.begin(), fractions.end(), double( 0 ),
336 std::cout <<
"PHG4MicromegasHitReco::process_event - sum: " <<
sum << std::endl;
343 for(
const auto& pair: fractions )
345 const int strip = pair.first;
346 if( strip < 0 || strip >= (
int) layergeom->get_strip_count( tileid ) )
continue;
348 const auto it = total_charges.lower_bound( strip );
349 if(
it != total_charges.end() &&
it->first == strip )
it->second += pair.second*gain;
350 else total_charges.insert(
it, std::make_pair( strip, pair.second*gain ) );
357 for(
const auto pair:total_charges )
360 const int strip = pair.first;
364 auto hit = hitset_it->second->getHit(hitkey);
369 hitset_it->second->addHitSpecificKey(hitkey, hit);
376 hittruthassoc->addAssoc(hitsetkey, hitkey, g4hit_it->first);
402 static constexpr
double Ar_dEdx = 2.44;
403 static constexpr
double iC4H10_dEdx = 5.93;
404 static constexpr
double mix_dEdx = 0.9*Ar_dEdx + 0.1*iC4H10_dEdx;
407 static constexpr
double Ar_ntot = 94;
408 static constexpr
double iC4H10_ntot = 195;
409 static constexpr
double mix_ntot = 0.9*Ar_ntot + 0.1*iC4H10_ntot;
412 static constexpr
double mix_electrons_per_gev = 1
e6*mix_ntot / mix_dEdx;
434 auto geonode_full = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename_full);
439 auto geonode_bare = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename_bare);
442 std::cout <<
PHWHERE <<
"Could not locate geometry node " << geonodename_bare << std::endl;
448 { std::cout <<
"PHG4MicromegasHitReco::setup_tiles - " <<
PHWHERE <<
" creating node " << geonodename_full << std::endl; }
454 runNode->addNode(newNode);
458 for(
auto iter = range.first; iter != range.second; ++iter )
460 const auto layer = iter->first;
468 for(
auto iter = range.first; iter != range.second; ++iter )
479 const bool is_first( iter == range.first );
480 cylinder->set_segmentation_type( is_first ?
489 cylinder->set_drift_direction( is_first ?
490 MicromegasDefs::DriftDirection::OUTWARD :
491 MicromegasDefs::DriftDirection::INWARD );
495 cylinder->set_pitch( is_first ? 25./256 : 50./256 );
498 { cylinder->identify( std::cout ); }
522 const TVector3& local_coords,
533 const auto pitch = layergeom->
get_pitch();
537 const int nstrips = 5.*sigma/pitch + 1;
538 const auto stripnum_min = std::clamp<int>( stripnum - nstrips, 0, strip_count );
539 const auto stripnum_max = std::clamp<int>( stripnum + nstrips, 0, strip_count );
545 for(
int strip = stripnum_min; strip <= stripnum_max; ++strip )
556 strip_location.x() - local_coords.x():
557 strip_location.z() - local_coords.z();
561 get_zigzag_fraction( xloc, sigma, pitch ):
562 get_rectangular_fraction( xloc, sigma, pitch );
565 charge_list.push_back( std::make_pair( strip, fraction ) );