ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4MicromegasHitReco.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4MicromegasHitReco.cc
1 
7 
10 
11 #include <trackbase/TrkrDefs.h>
12 #include <trackbase/TrkrHitv2.h>
13 #include <trackbase/TrkrHitSet.h>
16 
19 
20 #include <g4main/PHG4Hit.h>
21 #include <g4main/PHG4Hitv1.h>
23 
24 #include <phparameter/PHParameterInterface.h> // for PHParameterInterface
25 
27 #include <fun4all/SubsysReco.h> // for SubsysReco
28 
29 #include <phool/PHCompositeNode.h>
30 #include <phool/PHIODataNode.h>
31 #include <phool/PHNode.h>
32 #include <phool/PHNodeIterator.h>
33 #include <phool/PHRandomSeed.h>
34 #include <phool/getClass.h>
35 #include <phool/phool.h>
36 #include <phool/PHObject.h> // for PHObject
37 
38 
39 #include <TVector3.h>
40 
41 #include <gsl/gsl_randist.h>
42 #include <gsl/gsl_rng.h> // for gsl_rng_alloc
43 
44 #include <cassert>
45 #include <cmath> // for atan2, sqrt, M_PI
46 #include <cstdlib> // for exit
47 #include <iostream> // for operator<<, basic...
48 #include <map> // for _Rb_tree_const_it...
49 #include <numeric>
50 
51 namespace
52 {
54  template<class T>
55  inline constexpr T square( const T& x ) { return x*x; }
56 
58  template<class T>
59  inline T gaus( const T& x, const T& sigma )
60  { return std::exp( -square(x/sigma)/2 )/(sigma*std::sqrt(2*M_PI)); }
61 
63  template<class T>
64  inline T bind_angle( const T& angle )
65  {
66  if( angle >= M_PI ) return angle - 2*M_PI;
67  else if( angle < -M_PI ) return angle + 2*M_PI;
68  else return angle;
69  }
70 
71  // this corresponds to integrating a gaussian centered on zero and of width sigma from xloc - pitch/2 to xloc+pitch/2
72  template<class T>
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; }
75 
76  /*
77  this corresponds to integrating a gaussian centered on zero and of width sigma
78  convoluted with a zig-zag strip response function, which is triangular from xloc-pitch to xloc+pitch, with a maximum of 1 at xloc
79  */
80  template<class T>
81  inline T get_zigzag_fraction( const T& xloc, const T& sigma, const T& pitch )
82  {
83  return
84  // rising edge
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
87 
88  // descending edge
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;
91  }
92 
93 }
94 
95 //___________________________________________________________________________
96 PHG4MicromegasHitReco::PHG4MicromegasHitReco(const std::string &name, const std::string& detector)
97  : SubsysReco(name)
98  , PHParameterInterface(name)
99  , m_detector(detector)
100 {
101  // initialize rng
102  const uint seed = PHRandomSeed();
103  m_rng.reset( gsl_rng_alloc(gsl_rng_mt19937) );
104  gsl_rng_set( m_rng.get(), seed );
105 
107 }
108 
109 //___________________________________________________________________________
111 {
112 
114 
115  // load parameters
116  m_tmin = get_double_param("micromegas_tmin" );
117  m_tmax = get_double_param("micromegas_tmax" );
118  m_electrons_per_gev = get_double_param("micromegas_electrons_per_gev" );
119  m_gain = get_double_param("micromegas_gain");
120  m_cloud_sigma = get_double_param("micromegas_cloud_sigma");
121  m_diffusion_trans = get_double_param("micromegas_diffusion_trans");
122  m_zigzag_strips = get_int_param("micromegas_zigzag_strips");
123 
124  // printout
125  if( Verbosity() )
126  {
127  std::cout
128  << "PHG4MicromegasHitReco::InitRun\n"
129  << " m_tmin: " << m_tmin << "ns, m_tmax: " << m_tmax << "ns\n"
130  << " m_electrons_per_gev: " << m_electrons_per_gev << "\n"
131  << " m_gain: " << m_gain << "\n"
132  << " m_cloud_sigma: " << m_cloud_sigma << "cm\n"
133  << " m_diffusion_trans: " << m_diffusion_trans << "cm/sqrt(cm)\n"
134  << " m_zigzag_strips: " << std::boolalpha << m_zigzag_strips << "\n"
135  << std::endl;
136  }
137 
138  // setup tiles
139  setup_tiles( topNode );
140 
141  // get dst node
142  PHNodeIterator iter(topNode);
143  auto dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
144  if (!dstNode)
145  {
146  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
147  exit(1);
148  }
149 
150  // create hitset container if needed
151  auto hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
152  if (!hitsetcontainer)
153  {
154  std::cout << PHWHERE << "creating TRKR_HITSET." << std::endl;
155 
156  // find or create TRKR node
157  PHNodeIterator dstiter(dstNode);
158  auto trkrnode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
159  if (!trkrnode)
160  {
161  trkrnode = new PHCompositeNode("TRKR");
162  dstNode->addNode(trkrnode);
163  }
164 
165  // create container and add to the tree
166  hitsetcontainer = new TrkrHitSetContainerv1;
167  auto newNode = new PHIODataNode<PHObject>(hitsetcontainer, "TRKR_HITSET", "PHObject");
168  trkrnode->addNode(newNode);
169  }
170 
171  // create hit truth association if needed
172  auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
173  if (!hittruthassoc)
174  {
175  std::cout << PHWHERE << "creating TRKR_HITTRUTHASSOC." << std::endl;
176 
177  // find or create TRKR node
178  PHNodeIterator dstiter(dstNode);
179  auto trkrnode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
180  if (!trkrnode)
181  {
182  trkrnode = new PHCompositeNode("TRKR");
183  dstNode->addNode(trkrnode);
184  }
185 
186  hittruthassoc = new TrkrHitTruthAssocv1;
187  auto newNode = new PHIODataNode<PHObject>(hittruthassoc, "TRKR_HITTRUTHASSOC", "PHObject");
188  trkrnode->addNode(newNode);
189  }
190 
192 }
193 
194 //___________________________________________________________________________
196 {
197  // load relevant nodes
198  // G4Hits
199  const std::string g4hitnodename = "G4HIT_" + m_detector;
200  auto g4hitcontainer = findNode::getClass<PHG4HitContainer>(topNode, g4hitnodename);
201  assert(g4hitcontainer);
202 
203  // geometry
204  const auto geonodename = full_geonodename();
205  auto geonode = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename.c_str());
206  assert(geonode);
207 
208  // Get the TrkrHitSetContainer node
209  auto trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
210  assert(trkrhitsetcontainer);
211 
212  // Get the TrkrHitTruthAssoc node
213  auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
214  assert(hittruthassoc);
215 
216  // loop over layers in the g4hit container
217  auto layer_range = g4hitcontainer->getLayers();
218  for( auto layer_it = layer_range.first; layer_it != layer_range.second; ++layer_it )
219  {
220  // get layer
221  const auto layer = *layer_it;
222 
223  // get relevant geometry
224  auto layergeom = dynamic_cast<CylinderGeomMicromegas*>(geonode->GetLayerGeom(layer));
225  assert( layergeom );
226 
227  /*
228  * get the position of the detector mesh in local coordinates
229  * in local coordinate the mesh is a plane perpendicular to the y axis
230  * Its position along z depends on the drift direction
231  * it is used to calculate the drift distance of the primary electrons, and the
232  * corresponding transverse diffusion
233  */
234  const auto mesh_local_y = layergeom->get_drift_direction() == MicromegasDefs::DriftDirection::OUTWARD ?
235  layergeom->get_thickness()/2:
236  -layergeom->get_thickness()/2;
237 
238 // /*
239 // * get the radius of the detector mesh. It depends on the drift direction
240 // * it is used to calculate the drift distance of the primary electrons, and the
241 // * corresponding transverse diffusion
242 // */
243 // const auto mesh_radius = layergeom->get_drift_direction() == MicromegasDefs::DriftDirection::OUTWARD ?
244 // (layergeom->get_radius() + layergeom->get_thickness()/2):
245 // (layergeom->get_radius() - layergeom->get_thickness()/2);
246 
247  // get hits
248  const PHG4HitContainer::ConstRange g4hit_range = g4hitcontainer->getHits(layer);
249 
250  // loop over hits
251  for( auto g4hit_it = g4hit_range.first; g4hit_it != g4hit_range.second; ++g4hit_it )
252  {
253 
254  // get hit
255  PHG4Hit* g4hit = g4hit_it->second;
256 
257  // check time window
258  if(g4hit->get_t(0) > m_tmax) continue;
259  if(g4hit->get_t(1) < m_tmin) continue;
260 
261  // get world coordinates
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) );
264 
265  // make sure that the mid point is in one of the tiles
266  /*
267  * at this point we do not check the strip validity.
268  * This will be done when actually distributing electrons along the G4Hit track segment
269  */
270  const int tileid = layergeom->find_tile_cylindrical( (world_in+world_out)*0.5 );
271  if( tileid < 0 ) continue;
272 
273  /*
274  * in geant4 hits are generated on a cylinder located at layergeom->get_radius()
275  * the actual micromegas tiles however are plane surfaces, tengent to said cylinder
276  * one must convert the g4hit in and out positions from the cylinder back to the actual
277  * micromegas surface
278  */
279 
280  // make a local copy of the g4hit
281  PHG4Hitv1 g4hit_copy( g4hit );
282 
283  /*
284  * move hit coordinates to plane,
285  * update world coordinates
286  * get corresponding local coordinates
287  */
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) );
291 
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 );
294 
295  // number of primary elections
296  const auto nprimary = get_primary_electrons( &g4hit_copy );
297  if( !nprimary ) continue;
298 
299  // create hitset
300  const TrkrDefs::hitsetkey hitsetkey = MicromegasDefs::genHitSetKey( layer, layergeom->get_segmentation_type(), tileid );
301  const auto hitset_it = trkrhitsetcontainer->findOrAddHitSet(hitsetkey);
302 
303  // keep track of all charges
304  using charge_map_t = std::map<int,double>;
305  charge_map_t total_charges;
306 
307  // loop over primaries
308  for( uint ie = 0; ie < nprimary; ++ie )
309  {
310  // put the electron at a random position along the g4hit path
311  // in local reference frame, drift occurs along the y axis, from local y to mesh_local_y
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);
314 
315  if( m_diffusion_trans > 0 )
316  {
317  // add transeverse diffusion
318  // first convert to polar coordinates
319  const double y = local.y();
320  const double drift_distance = std::abs(y - mesh_local_y);
321  const double diffusion = gsl_ran_gaussian(m_rng.get(), m_diffusion_trans*std::sqrt(drift_distance));
322  const double diffusion_angle = gsl_ran_flat(m_rng.get(), -M_PI, M_PI);
323 
324  // diffusion occurs in x,z plane with a magnitude 'diffusion' and an angle 'diffusion angle'
325  local += TVector3( diffusion*std::cos(diffusion_angle), 0, diffusion*std::sin(diffusion_angle) );
326  }
327 
328  // distribute charge among adjacent strips
329  const auto fractions = distribute_charge( layergeom, tileid, local, m_cloud_sigma );
330 
331  // make sure fractions adds up to unity
332  if( Verbosity() > 0 )
333  {
334  const auto sum = std::accumulate( fractions.begin(), fractions.end(), double( 0 ),
335  []( double value, const charge_pair_t& pair ) { return value + pair.second; } );
336  std::cout << "PHG4MicromegasHitReco::process_event - sum: " << sum << std::endl;
337  }
338 
339  // generate gain for this electron
340  const auto gain = get_single_electron_amplification();
341 
342  // merge to total charges
343  for( const auto& pair: fractions )
344  {
345  const int strip = pair.first;
346  if( strip < 0 || strip >= (int) layergeom->get_strip_count( tileid ) ) continue;
347 
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 ) );
351  }
352 
353  }
354 
355  // generate the key for this hit
356  // loop over strips in list
357  for( const auto pair:total_charges )
358  {
359  // get strip and bound check
360  const int strip = pair.first;
361 
362  // get hit from hitset
364  auto hit = hitset_it->second->getHit(hitkey);
365  if( !hit )
366  {
367  // create hit and insert in hitset
368  hit = new TrkrHitv2;
369  hitset_it->second->addHitSpecificKey(hitkey, hit);
370  }
371 
372  // add energy from g4hit
373  hit->addEnergy( pair.second );
374 
375  // associate this hitset and hit to the geant4 hit key
376  hittruthassoc->addAssoc(hitsetkey, hitkey, g4hit_it->first);
377  }
378 
379  }
380 
381  }
382 
384 }
385 
386 //___________________________________________________________________________
388 {
389  // default timing window (ns)
390  /*
391  * see https://indico.bnl.gov/event/8548/contributions/37753/attachments/28212/43343/2020_05_Proposal_sPhenixMonitoring_update_19052020.pptx slide 10
392  * small negative time for tmin is set to catch out of time, same-bunch pileup events
393  * similar value is used in PHG4InttReco
394  */
395  set_default_double_param("micromegas_tmin", -20 );
396  set_default_double_param("micromegas_tmax", 800 );
397 
398  // gas data from
399  // http://www.slac.stanford.edu/pubs/icfa/summer98/paper3/paper3.pdf
400  // assume Ar/iC4H10 90/10, at 20C and 1atm
401  // dedx (KeV/cm) for MIP
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;
405 
406  // number of electrons per MIP (cm-1)
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;
410 
411  // number of electrons per gev
412  static constexpr double mix_electrons_per_gev = 1e6*mix_ntot / mix_dEdx;
413  set_default_double_param("micromegas_electrons_per_gev", mix_electrons_per_gev );
414 
415  // gain
416  set_default_double_param("micromegas_gain", 2000 );
417 
418  // electron cloud sigma, after avalanche (cm)
419  set_default_double_param("micromegas_cloud_sigma", 0.04 );
420 
421  // transverse diffusion (cm/sqrt(cm))
422  set_default_double_param("micromegas_diffusion_trans", 0.03 );
423 
424  // zigzag strips
425  set_default_int_param("micromegas_zigzag_strips", true );
426 }
427 
428 //___________________________________________________________________________
430 {
431 
432  // get geometry
433  const auto geonodename_full = full_geonodename();
434  auto geonode_full = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename_full);
435  if (!geonode_full)
436  {
437  // if full geometry (cylinder + tiles) do not exist, try create it from bare geometry (cylinder only)
438  const auto geonodename_bare = bare_geonodename();
439  auto geonode_bare = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename_bare);
440  if( !geonode_bare )
441  {
442  std::cout << PHWHERE << "Could not locate geometry node " << geonodename_bare << std::endl;
443  exit(1);
444  }
445 
446  // create new node
447  if( Verbosity() )
448  { std::cout << "PHG4MicromegasHitReco::setup_tiles - " << PHWHERE << " creating node " << geonodename_full << std::endl; }
449 
450  geonode_full = new PHG4CylinderGeomContainer();
451  PHNodeIterator iter(topNode);
452  auto runNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "RUN"));
453  auto newNode = new PHIODataNode<PHObject>(geonode_full, geonodename_full, "PHObject");
454  runNode->addNode(newNode);
455 
456  // copy cylinders
457  PHG4CylinderGeomContainer::ConstRange range = geonode_bare->get_begin_end();
458  for( auto iter = range.first; iter != range.second; ++iter )
459  {
460  const auto layer = iter->first;
461  const auto cylinder = static_cast<CylinderGeomMicromegas*>(iter->second);
462  geonode_full->AddLayerGeom( layer, new CylinderGeomMicromegas( *cylinder ) );
463  }
464  }
465 
466  // get cylinders
467  PHG4CylinderGeomContainer::ConstRange range = geonode_full->get_begin_end();
468  for( auto iter = range.first; iter != range.second; ++iter )
469  {
470  auto cylinder = static_cast<CylinderGeomMicromegas*>(iter->second);
471 
472  // assign tiles
473  cylinder->set_tiles( m_tiles );
474 
475  /*
476  * asign segmentation type and pitch
477  * assume first layer in phi, other(s) are z
478  */
479  const bool is_first( iter == range.first );
480  cylinder->set_segmentation_type( is_first ?
483 
484  /*
485  * assign drift direction
486  * assume first layer is outward, with readout plane at the top, and others are inward, with readout plane at the bottom
487  * this is used to properly implement transverse diffusion in ::process_event
488  */
489  cylinder->set_drift_direction( is_first ?
490  MicromegasDefs::DriftDirection::OUTWARD :
491  MicromegasDefs::DriftDirection::INWARD );
492 
493  // pitch
494  /* they correspond to 256 channels along the phi direction, and 256 along the z direction, assuming 25x50 tiles */
495  cylinder->set_pitch( is_first ? 25./256 : 50./256 );
496 
497  if( Verbosity() )
498  { cylinder->identify( std::cout ); }
499 
500  }
501 }
502 
503 //___________________________________________________________________________
505 { return gsl_ran_poisson(m_rng.get(), g4hit->get_eion() * m_electrons_per_gev); }
506 
507 //___________________________________________________________________________
509 {
510  /*
511  * to handle gain fluctuations, an exponential distribution is used, similar to what used for the GEMS
512  * (simulations/g4simulations/g4tpc/PHG4TpcPadPlaneReadout::getSingleEGEMAmplification)
513  * One must get a different random number for each primary electron for this to be valid
514  */
515  return gsl_ran_exponential(m_rng.get(), m_gain);
516 }
517 
518 //___________________________________________________________________________
520  CylinderGeomMicromegas* layergeom,
521  uint tileid,
522  const TVector3& local_coords,
523  double sigma ) const
524 {
525 
526  // find tile and strip matching center position
527  auto stripnum = layergeom->find_strip_from_local_coords( tileid, local_coords );
528 
529  // check tile and strip
530  if( stripnum < 0 ) return charge_list_t();
531 
532  // store pitch and radius
533  const auto pitch = layergeom->get_pitch();
534 
535  // find relevant strip indices
536  const auto strip_count = layergeom->get_strip_count( tileid );
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 );
540 
541  // prepare charge list
542  charge_list_t charge_list;
543 
544  // loop over strips
545  for( int strip = stripnum_min; strip <= stripnum_max; ++strip )
546  {
547  // get strip center
548  const TVector3 strip_location = layergeom->get_local_coordinates( tileid, strip );
549 
550  /*
551  * find relevant strip coordinate with respect to location
552  * in local coordinate, phi segmented view has strips along z and measures along x
553  * in local coordinate, z segmented view has strips along phi and measures along z
554  */
556  strip_location.x() - local_coords.x():
557  strip_location.z() - local_coords.z();
558 
559  // calculate charge fraction
560  const auto fraction = m_zigzag_strips ?
561  get_zigzag_fraction( xloc, sigma, pitch ):
562  get_rectangular_fraction( xloc, sigma, pitch );
563 
564  // store
565  charge_list.push_back( std::make_pair( strip, fraction ) );
566  }
567 
568  return charge_list;
569 }