ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4TpcElectronDrift.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4TpcElectronDrift.cc
1 // this is the new containers version
2 // it uses the same MapToPadPlane as the old containers version
3 
4 #include "PHG4TpcElectronDrift.h"
5 #include "PHG4TpcDistortion.h"
6 #include "PHG4TpcPadPlane.h" // for PHG4TpcPadPlane
7 
8 #include <g4main/PHG4Hit.h>
10 
11 #include <trackbase/TrkrDefs.h>
12 #include <trackbase/TrkrHit.h> // for TrkrHit
13 #include <trackbase/TrkrHitSet.h>
15 #include <trackbase/TrkrHitTruthAssoc.h> // for TrkrHitTruthA...
17 #include <trackbase/TrkrHitv2.h>
18 
19 #include <tpc/TpcDefs.h>
20 
22 
23 #include <phparameter/PHParameterInterface.h> // for PHParameterIn...
24 #include <phparameter/PHParameters.h>
25 #include <phparameter/PHParametersContainer.h>
26 
27 #include <pdbcalbase/PdbParameterMapContainer.h>
28 
30 #include <fun4all/Fun4AllServer.h>
31 #include <fun4all/SubsysReco.h> // for SubsysReco
32 
33 #include <phool/PHCompositeNode.h>
34 #include <phool/PHDataNode.h> // for PHDataNode
35 #include <phool/PHIODataNode.h>
36 #include <phool/PHNode.h> // for PHNode
37 #include <phool/PHNodeIterator.h>
38 #include <phool/PHObject.h> // for PHObject
39 #include <phool/PHRandomSeed.h>
40 #include <phool/getClass.h>
41 #include <phool/phool.h> // for PHWHERE
42 
43 #include <TFile.h>
44 #include <TH1.h>
45 #include <TH2.h>
46 #include <TNtuple.h>
47 #include <TSystem.h>
48 
49 #include <gsl/gsl_randist.h>
50 #include <gsl/gsl_rng.h> // for gsl_rng_alloc
51 
52 #include <cassert>
53 #include <cmath> // for sqrt, abs, NAN
54 #include <cstdlib> // for exit
55 #include <iostream>
56 #include <map> // for _Rb_tree_cons...
57 #include <utility> // for pair
58 
59 namespace
60 {
61  template <class T>
62  inline constexpr T square(const T &x)
63  {
64  return x * x;
65  }
66 } // namespace
67 
69  : SubsysReco(name)
70  , PHParameterInterface(name)
71  , temp_hitsetcontainer(new TrkrHitSetContainerv1)
72  , single_hitsetcontainer(new TrkrHitSetContainerv1)
73 {
75  RandomGenerator.reset(gsl_rng_alloc(gsl_rng_mt19937));
77 }
78 
79 //_____________________________________________________________
81 {
82  padplane->Init(topNode);
83  event_num = 0;
85 }
86 
87 //_____________________________________________________________
89 {
90  PHNodeIterator iter(topNode);
91 
92  // Looking for the DST node
93  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
94  if (!dstNode)
95  {
96  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
97  exit(1);
98  }
99  auto runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
100  auto parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PAR"));
101  const std::string paramnodename = "G4CELLPARAM_" + detector;
102  const std::string geonodename = "G4CELLPAR_" + detector;
103  const std::string tpcgeonodename = "G4GEO_" + detector;
104  hitnodename = "G4HIT_" + detector;
105  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
106  if (!g4hit)
107  {
108  std::cout << Name() << " Could not locate G4HIT node " << hitnodename << std::endl;
109  topNode->print();
110  gSystem->Exit(1);
111  exit(1);
112  }
113 
114  // new containers
115  hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
116  if (!hitsetcontainer)
117  {
118  PHNodeIterator dstiter(dstNode);
119  auto DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
120  if (!DetNode)
121  {
122  DetNode = new PHCompositeNode("TRKR");
123  dstNode->addNode(DetNode);
124  }
125 
127  auto newNode = new PHIODataNode<PHObject>(hitsetcontainer, "TRKR_HITSET", "PHObject");
128  DetNode->addNode(newNode);
129  }
130 
131  hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
132  if (!hittruthassoc)
133  {
134  PHNodeIterator dstiter(dstNode);
135  auto DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
136  if (!DetNode)
137  {
138  DetNode = new PHCompositeNode("TRKR");
139  dstNode->addNode(DetNode);
140  }
141 
143  auto newNode = new PHIODataNode<PHObject>(hittruthassoc, "TRKR_HITTRUTHASSOC", "PHObject");
144  DetNode->addNode(newNode);
145  }
146 
147  seggeonodename = "CYLINDERCELLGEOM_SVTX"; // + detector;
148  PHG4CylinderCellGeomContainer *seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, seggeonodename.c_str());
149  if (!seggeo)
150  {
151  seggeo = new PHG4CylinderCellGeomContainer();
152  auto runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
153  auto newNode = new PHIODataNode<PHObject>(seggeo, seggeonodename.c_str(), "PHObject");
154  runNode->addNode(newNode);
155  }
156 
158  PHNodeIterator runIter(runNode);
159  auto RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", detector));
160  if (!RunDetNode)
161  {
162  RunDetNode = new PHCompositeNode(detector);
163  runNode->addNode(RunDetNode);
164  }
165  SaveToNodeTree(RunDetNode, paramnodename);
166 
167  // save this to the parNode for use
168  PHNodeIterator parIter(parNode);
169  auto ParDetNode = dynamic_cast<PHCompositeNode *>(parIter.findFirst("PHCompositeNode", detector));
170  if (!ParDetNode)
171  {
172  ParDetNode = new PHCompositeNode(detector);
173  parNode->addNode(ParDetNode);
174  }
175  PutOnParNode(ParDetNode, geonodename);
176 
177  // find Tpc Geo
178  PHNodeIterator tpcpariter(ParDetNode);
179  auto tpcparams = findNode::getClass<PHParametersContainer>(ParDetNode, tpcgeonodename);
180  if (!tpcparams)
181  {
182  const std::string runparamname = "G4GEOPARAM_" + detector;
183  auto tpcpdbparams = findNode::getClass<PdbParameterMapContainer>(RunDetNode, runparamname);
184  if (tpcpdbparams)
185  {
186  tpcparams = new PHParametersContainer(detector);
187  if (Verbosity()) tpcpdbparams->print();
188  tpcparams->CreateAndFillFrom(tpcpdbparams, detector);
189  ParDetNode->addNode(new PHDataNode<PHParametersContainer>(tpcparams, tpcgeonodename));
190  }
191  else
192  {
193  std::cout << "PHG4TpcElectronDrift::InitRun - failed to find " << runparamname << " in order to initialize " << tpcgeonodename << ". Aborting run ..." << std::endl;
195  }
196  }
197  assert(tpcparams);
198 
199  if (Verbosity()) tpcparams->Print();
200  const PHParameters *tpcparam = tpcparams->GetParameters(0);
201  assert(tpcparam);
202  tpc_length = tpcparam->get_double_param("tpc_length");
203 
204  diffusion_long = get_double_param("diffusion_long");
205  added_smear_sigma_long = get_double_param("added_smear_long");
206  diffusion_trans = get_double_param("diffusion_trans");
207  added_smear_sigma_trans = get_double_param("added_smear_trans");
208  drift_velocity = get_double_param("drift_velocity");
209  min_time = 0.0;
210  max_time = (tpc_length / 1.75) / drift_velocity;
211  electrons_per_gev = get_double_param("electrons_per_gev");
212  min_active_radius = get_double_param("min_active_radius");
213  max_active_radius = get_double_param("max_active_radius");
214 
215  auto se = Fun4AllServer::instance();
216  dlong = new TH1F("difflong", "longitudinal diffusion", 100, diffusion_long - diffusion_long / 2., diffusion_long + diffusion_long / 2.);
217  se->registerHisto(dlong);
218  dtrans = new TH1F("difftrans", "transversal diffusion", 100, diffusion_trans - diffusion_trans / 2., diffusion_trans + diffusion_trans / 2.);
219  se->registerHisto(dtrans);
220 
221  do_ElectronDriftQAHistos = false; // Whether or not to produce an ElectronDriftQA.root file with useful info
223  {
224  hitmapstart = new TH2F("hitmapstart", "g4hit starting X-Y locations", 1560, -78, 78, 1560, -78, 78);
225  hitmapend = new TH2F("hitmapend", "g4hit final X-Y locations", 1560, -78, 78, 1560, -78, 78);
226  z_startmap = new TH2F("z_startmap", "g4hit starting Z vs. R locations", 2000, -100, 100, 780, 0, 78);
227  deltaphi = new TH2F("deltaphi", "Total delta phi; phi (rad);#Delta phi (rad)", 600, -M_PI, M_PI, 1000, -.2, .2);
228  deltaRphinodiff = new TH2F("deltaRphinodiff", "Total delta R*phi, no diffusion; r (cm);#Delta R*phi (cm)", 600, 20, 80, 1000, -3, 5);
229  deltaphivsRnodiff = new TH2F("deltaphivsRnodiff", "Total delta phi vs. R; phi (rad);#Delta phi (rad)", 600, 20, 80, 1000, -.2, .2);
230  deltaz = new TH2F("deltaz", "Total delta z; z (cm);#Delta z (cm)", 1000, 0, 100, 1000, -.5, 5);
231  deltaphinodiff = new TH2F("deltaphinodiff", "Total delta phi (no diffusion, only SC distortion); phi (rad);#Delta phi (rad)", 600, -M_PI, M_PI, 1000, -.2, .2);
232  deltaphinodist = new TH2F("deltaphinodist", "Total delta phi (no SC distortion, only diffusion); phi (rad);#Delta phi (rad)", 600, -M_PI, M_PI, 1000, -.2, .2);
233  deltar = new TH2F("deltar", "Total Delta r; r (cm);#Delta r (cm)", 580, 20, 78, 1000, -3, 5);
234  deltarnodiff = new TH2F("deltarnodiff", "Delta r (no diffusion, only SC distortion); r (cm);#Delta r (cm)", 580, 20, 78, 1000, -2, 5);
235  deltarnodist = new TH2F("deltarnodist", "Delta r (no SC distortion, only diffusion); r (cm);#Delta r (cm)", 580, 20, 78, 1000, -2, 5);
236  }
237 
238  if (Verbosity())
239  {
240  // eval tree only when verbosity is on
241  m_outf.reset(new TFile("nt_out.root", "recreate"));
242  nt = new TNtuple("nt", "electron drift stuff", "hit:ts:tb:tsig:rad:zstart:zfinal");
243  nthit = new TNtuple("nthit", "TrkrHit collecting", "layer:phipad:zbin:neffelectrons");
244  ntfinalhit = new TNtuple("ntfinalhit", "TrkrHit collecting", "layer:phipad:zbin:neffelectrons");
245  ntpad = new TNtuple("ntpad", "electron by electron pad centroid", "layer:phigem:phiclus:zgem:zclus");
246  se->registerHisto(nt);
247  se->registerHisto(nthit);
248  se->registerHisto(ntpad);
249  }
250  padplane->InitRun(topNode);
251  padplane->CreateReadoutGeometry(topNode, seggeo);
252 
254 }
255 
257 {
258  static constexpr unsigned int print_layer = 18;
259 
260  // tells m_distortionMap which event to look at
261  if (m_distortionMap)
262  {
263  m_distortionMap->load_event(event_num);
264  }
265 
266  // g4hits
267  auto g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
268  if (!g4hit)
269  {
270  std::cout << "Could not locate g4 hit node " << hitnodename << std::endl;
271  gSystem->Exit(1);
272  }
273 
274  PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits();
275  //std::cout << "g4hits size " << g4hit->size() << std::endl;
276  unsigned int count_g4hits = 0;
277  int count_electrons = 0;
278 
279  double ecollectedhits = 0.0;
280  int ncollectedhits = 0;
281  double ihit = 0;
282  unsigned int dump_interval = 5000; // dump temp_hitsetcontainer to the node tree after this many g4hits
283  unsigned int dump_counter = 0;
284  for (auto hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
285  {
286  count_g4hits++;
287  dump_counter++;
288 
289  const double t0 = fmax(hiter->second->get_t(0), hiter->second->get_t(1));
290  if (t0 > max_time)
291  {
292  continue;
293  }
294 
295  // for very high occupancy events, accessing the TrkrHitsets on the node tree for every drifted electron seems to be very slow
296  // Instead, use a temporary map to accumulate the charge from all drifted electrons, then copy to the node tree later
297  double eion = hiter->second->get_eion();
298  unsigned int n_electrons = gsl_ran_poisson(RandomGenerator.get(), eion * electrons_per_gev);
299  count_electrons += n_electrons;
300 
301  /*
302  if(count_g4hits%50000 == 0)
303  std::cout << " g4hit->size() " << g4hit->size() << " count_g4hits " << count_g4hits << " remaining " <<
304  g4hit->size() - count_g4hits << " count_electrons " << count_electrons << std::endl;
305  */
306 
307  if (Verbosity() > 100)
308  std::cout << " new hit with t0, " << t0 << " g4hitid " << hiter->first
309  << " eion " << eion << " n_electrons " << n_electrons
310  << " entry z " << hiter->second->get_z(0) << " exit z " << hiter->second->get_z(1) << " avg z" << (hiter->second->get_z(0) + hiter->second->get_z(1)) / 2.0
311  << std::endl;
312 
313  if (n_electrons == 0)
314  {
315  continue;
316  }
317 
318  if (Verbosity() > 100)
319  {
320  std::cout << std::endl
321  << "electron drift: g4hit " << hiter->first << " created electrons: " << n_electrons
322  << " from " << eion * 1000000 << " keV" << std::endl;
323  std::cout << " entry x,y,z = " << hiter->second->get_x(0) << " " << hiter->second->get_y(0) << " " << hiter->second->get_z(0)
324  << " radius " << sqrt(pow(hiter->second->get_x(0), 2) + pow(hiter->second->get_y(0), 2)) << std::endl;
325  std::cout << " exit x,y,z = " << hiter->second->get_x(1) << " " << hiter->second->get_y(1) << " " << hiter->second->get_z(1)
326  << " radius " << sqrt(pow(hiter->second->get_x(1), 2) + pow(hiter->second->get_y(1), 2)) << std::endl;
327  }
328 
329  for (unsigned int i = 0; i < n_electrons; i++)
330  {
331  // We choose the electron starting position at random from a flat distribution along the path length
332  // the parameter t is the fraction of the distance along the path betwen entry and exit points, it has values between 0 and 1
333  const double f = gsl_ran_flat(RandomGenerator.get(), 0.0, 1.0);
334 
335  const double x_start = hiter->second->get_x(0) + f * (hiter->second->get_x(1) - hiter->second->get_x(0));
336  const double y_start = hiter->second->get_y(0) + f * (hiter->second->get_y(1) - hiter->second->get_y(0));
337  const double z_start = hiter->second->get_z(0) + f * (hiter->second->get_z(1) - hiter->second->get_z(0));
338  const double t_start = hiter->second->get_t(0) + f * (hiter->second->get_t(1) - hiter->second->get_t(0));
339 
340  const double r_sigma = diffusion_trans * sqrt(tpc_length / 2. - std::abs(z_start));
341  const double rantrans =
342  gsl_ran_gaussian(RandomGenerator.get(), r_sigma) +
343  gsl_ran_gaussian(RandomGenerator.get(), added_smear_sigma_trans);
344 
345  const double t_path = (tpc_length / 2. - std::abs(z_start)) / drift_velocity;
346  const double t_sigma = diffusion_long * sqrt(tpc_length / 2. - std::abs(z_start)) / drift_velocity;
347  const double rantime =
348  gsl_ran_gaussian(RandomGenerator.get(), t_sigma) +
349  gsl_ran_gaussian(RandomGenerator.get(), added_smear_sigma_long) / drift_velocity;
350  const double t_final = t_start + t_path + rantime;
351  if (t_final < min_time || t_final > max_time) continue;
352 
353  double z_final;
354  if (z_start < 0)
355  z_final = -tpc_length / 2. + t_final * drift_velocity;
356  else
357  z_final = tpc_length / 2. - t_final * drift_velocity;
358 
359  const double radstart = std::sqrt(square(x_start) + square(y_start));
360  const double phistart = std::atan2(y_start, x_start);
361  const double ranphi = gsl_ran_flat(RandomGenerator.get(), -M_PI, M_PI);
362 
363  double x_final = x_start + rantrans * std::cos(ranphi); // Initialize these to be only diffused first, will be overwritten if doing SC distortion
364  double y_final = y_start + rantrans * std::sin(ranphi);
365 
366  double rad_final = sqrt(square(x_final) + square(y_final));
367  double phi_final = atan2(y_final, x_final);
368 
370  {
371  z_startmap->Fill(z_start, radstart); // map of starting location in Z vs. R
372  deltaphinodist->Fill(phistart, rantrans / rad_final); // delta phi no distortion, just diffusion+smear
373  deltarnodist->Fill(radstart, rantrans); // delta r no distortion, just diffusion+smear
374  }
375 
376  if (m_distortionMap)
377  {
378  const double r_distortion = m_distortionMap->get_r_distortion(radstart, phistart, z_start);
379  const double phi_distortion = m_distortionMap->get_rphi_distortion(radstart, phistart, z_start)/radstart;
380  const double z_distortion = m_distortionMap->get_z_distortion(radstart, phistart, z_start);
381 
382  rad_final+=r_distortion;
383  phi_final+=phi_distortion;
384  z_final += z_distortion;
385 
386  x_final=rad_final*std::cos(phi_final);
387  y_final=rad_final*std::sin(phi_final);
388 
390  {
391  const double phi_final_nodiff = phistart+phi_distortion;
392  const double rad_final_nodiff = radstart+r_distortion;
393  deltarnodiff->Fill(radstart, rad_final_nodiff - radstart); //delta r no diffusion, just distortion
394  deltaphinodiff->Fill(phistart, phi_final_nodiff - phistart); //delta phi no diffusion, just distortion
395  deltaphivsRnodiff->Fill(radstart, phi_final_nodiff - phistart);
396  deltaRphinodiff->Fill(radstart, rad_final_nodiff * phi_final_nodiff - radstart * phistart);
397 
398  // Fill Diagnostic plots, written into ElectronDriftQA.root
399  hitmapstart->Fill(x_start, y_start); // G4Hit starting positions
400  hitmapend->Fill(x_final, y_final); //INcludes diffusion and distortion
401  deltar->Fill(radstart, rad_final - radstart); //total delta r
402  deltaphi->Fill(phistart, phi_final - phistart); // total delta phi
403  deltaz->Fill(z_start, z_distortion); // map of distortion in Z (time)
404  }
405  }
406 
407  // remove electrons outside of our acceptance. Careful though, electrons from just inside 30 cm can contribute in the 1st active layer readout, so leave a little margin
408  if (rad_final < min_active_radius - 2.0 || rad_final > max_active_radius + 1.0)
409  {
410  continue;
411  }
412 
413  if (Verbosity() > 1000)
414  {
415  std::cout << "electron " << i << " g4hitid " << hiter->first << " f " << f << std::endl;
416  std::cout << "radstart " << radstart << " x_start: " << x_start
417  << ", y_start: " << y_start
418  << ",z_start: " << z_start
419  << " t_start " << t_start
420  << " t_path " << t_path
421  << " t_sigma " << t_sigma
422  << " rantime " << rantime
423  << std::endl;
424 
425  //if( sqrt(x_start*x_start+y_start*y_start) > 68.0 && sqrt(x_start*x_start+y_start*y_start) < 72.0)
426  std::cout << " rad_final " << rad_final << " x_final " << x_final << " y_final " << y_final
427  << " z_final " << z_final << " t_final " << t_final << " zdiff " << z_final - z_start << std::endl;
428  }
429 
430  if (Verbosity() > 0)
431  {
432  assert(nt);
433  nt->Fill(ihit, t_start, t_final, t_sigma, rad_final, z_start, z_final);
434  }
435  // this fills the cells and updates the hits in temp_hitsetcontainer for this drifted electron hitting the GEM stack
436  MapToPadPlane(x_final, y_final, z_final, hiter, ntpad, nthit);
437  } // end loop over electrons for this g4hit
438 
439  // The hit-truth association has to be done for each g4hit
440  // we use the single_hitsetcontainer for this
441 
443  for (TrkrHitSetContainer::ConstIterator single_hitset_iter = single_hitset_range.first;
444  single_hitset_iter != single_hitset_range.second;
445  ++single_hitset_iter)
446  {
447  // we have an itrator to one TrkrHitSet for the Tpc from the single_hitsetcontainer
448  TrkrDefs::hitsetkey node_hitsetkey = single_hitset_iter->first;
449  const unsigned int layer = TrkrDefs::getLayer(node_hitsetkey);
450  const int sector = TpcDefs::getSectorId(node_hitsetkey);
451  const int side = TpcDefs::getSide(node_hitsetkey);
452 
453  if (Verbosity() > 2)
454  std::cout << " hitsetkey " << node_hitsetkey << " layer " << layer << " sector " << sector << " side " << side << std::endl;
455  // get all of the hits from the single hitset
456  TrkrHitSet::ConstRange single_hit_range = single_hitset_iter->second->getHits();
457  for (TrkrHitSet::ConstIterator single_hit_iter = single_hit_range.first;
458  single_hit_iter != single_hit_range.second;
459  ++single_hit_iter)
460  {
461  TrkrDefs::hitkey single_hitkey = single_hit_iter->first;
462 
463  // Add the hit-g4hit association
464  // no need to check for duplicates, since the hit is new
465  hittruthassoc->addAssoc(node_hitsetkey, single_hitkey, hiter->first);
466  if (Verbosity() > 100)
467  std::cout << " adding assoc for node_hitsetkey " << node_hitsetkey << " single_hitkey " << single_hitkey << " g4hitkey " << hiter->first << std::endl;
468  }
469  }
470 
471  // Dump the temp_hitsetcontainer to the node tree and reset it
472  // - after every "dump_interval" g4hits
473  // - if this is the last g4hit
474  if (dump_counter >= dump_interval || count_g4hits == g4hit->size())
475  {
476  //std::cout << " dump_counter " << dump_counter << " count_g4hits " << count_g4hits << std::endl;
477 
478  double eg4hit = 0.0;
480  for (TrkrHitSetContainer::ConstIterator temp_hitset_iter = temp_hitset_range.first;
481  temp_hitset_iter != temp_hitset_range.second;
482  ++temp_hitset_iter)
483  {
484  // we have an itrator to one TrkrHitSet for the Tpc from the temp_hitsetcontainer
485  TrkrDefs::hitsetkey node_hitsetkey = temp_hitset_iter->first;
486  const unsigned int layer = TrkrDefs::getLayer(node_hitsetkey);
487  const int sector = TpcDefs::getSectorId(node_hitsetkey);
488  const int side = TpcDefs::getSide(node_hitsetkey);
489  if (Verbosity() > 100)
490  std::cout << "PHG4TpcElectronDrift: temp_hitset with key: " << node_hitsetkey << " in layer " << layer
491  << " with sector " << sector << " side " << side << std::endl;
492 
493  // find or add this hitset on the node tree
494  TrkrHitSetContainer::Iterator node_hitsetit = hitsetcontainer->findOrAddHitSet(node_hitsetkey);
495 
496  // get all of the hits from the temporary hitset
497  TrkrHitSet::ConstRange temp_hit_range = temp_hitset_iter->second->getHits();
498  for (TrkrHitSet::ConstIterator temp_hit_iter = temp_hit_range.first;
499  temp_hit_iter != temp_hit_range.second;
500  ++temp_hit_iter)
501  {
502  TrkrDefs::hitkey temp_hitkey = temp_hit_iter->first;
503  TrkrHit *temp_tpchit = temp_hit_iter->second;
504  if (Verbosity() > 10 && layer == print_layer)
505  {
506  std::cout << " temp_hitkey " << temp_hitkey << " l;ayer " << layer << " pad " << TpcDefs::getPad(temp_hitkey)
507  << " z bin " << TpcDefs::getTBin(temp_hitkey)
508  << " energy " << temp_tpchit->getEnergy() << " eg4hit " << eg4hit << std::endl;
509 
510  eg4hit += temp_tpchit->getEnergy();
511  ecollectedhits += temp_tpchit->getEnergy();
512  ncollectedhits++;
513  }
514 
515  // find or add this hit to the node tree
516  TrkrHit *node_hit = node_hitsetit->second->getHit(temp_hitkey);
517  if (!node_hit)
518  {
519  // Otherwise, create a new one
520  node_hit = new TrkrHitv2();
521  node_hitsetit->second->addHitSpecificKey(temp_hitkey, node_hit);
522  }
523 
524  // Either way, add the energy to it
525  node_hit->addEnergy(temp_tpchit->getEnergy());
526 
527  } // end loop over temp hits
528 
529  if (Verbosity() > 100 && layer == print_layer)
530  std::cout << " ihit " << ihit << " collected energy = " << eg4hit << std::endl;
531 
532  } // end loop over temp hitsets
533 
534  // erase all entries in the temp hitsetcontainer
535  temp_hitsetcontainer->Reset();
536 
537  // reset the dump counter
538  dump_counter = 0;
539  } // end copy of temp hitsetcontainer to node tree hitsetcontainer
540 
541  ++ihit;
542 
543  // we are done with this until the next g4hit
544  single_hitsetcontainer->Reset();
545 
546  } // end loop over g4hits
547 
548  if (Verbosity() > 2)
549  {
550  std::cout << "From PHG4TpcElectronDrift: hitsetcontainer printout at end:" << std::endl;
551  // We want all hitsets for the Tpc
553  for (TrkrHitSetContainer::ConstIterator hitset_iter = hitset_range.first;
554  hitset_iter != hitset_range.second;
555  ++hitset_iter)
556  {
557  // we have an itrator to one TrkrHitSet for the Tpc from the trkrHitSetContainer
558  TrkrDefs::hitsetkey hitsetkey = hitset_iter->first;
559  const unsigned int layer = TrkrDefs::getLayer(hitsetkey);
560  if (layer != print_layer) continue;
561  const int sector = TpcDefs::getSectorId(hitsetkey);
562  const int side = TpcDefs::getSide(hitsetkey);
563 
564  std::cout << "PHG4TpcElectronDrift: hitset with key: " << hitsetkey << " in layer " << layer << " with sector " << sector << " side " << side << std::endl;
565 
566  // get all of the hits from this hitset
567  TrkrHitSet *hitset = hitset_iter->second;
568  TrkrHitSet::ConstRange hit_range = hitset->getHits();
569  for (TrkrHitSet::ConstIterator hit_iter = hit_range.first;
570  hit_iter != hit_range.second;
571  ++hit_iter)
572  {
573  TrkrDefs::hitkey hitkey = hit_iter->first;
574  TrkrHit *tpchit = hit_iter->second;
575  std::cout << " hitkey " << hitkey << " pad " << TpcDefs::getPad(hitkey) << " z bin " << TpcDefs::getTBin(hitkey)
576  << " energy " << tpchit->getEnergy() << std::endl;
577  }
578  }
579  }
580 
581  if (Verbosity() > 1000)
582  {
583  std::cout << "From PHG4TpcElectronDrift: hittruthassoc dump:" << std::endl;
585  }
586 
587  ++event_num; // if doing more than one event, event_num will be incremented.
589 }
590 
591 void PHG4TpcElectronDrift::MapToPadPlane(const double x_gem, const double y_gem, const double t_gem, PHG4HitContainer::ConstIterator hiter, TNtuple *ntpad, TNtuple *nthit)
592 {
593  padplane->MapToPadPlane(single_hitsetcontainer.get(), temp_hitsetcontainer.get(), hittruthassoc, x_gem, y_gem, t_gem, hiter, ntpad, nthit);
594 }
595 
597 {
598  if (Verbosity() > 0)
599  {
600  assert(m_outf);
601  assert(nt);
602  assert(ntpad);
603  assert(nthit);
604  assert(ntfinalhit);
605 
606  m_outf->cd();
607  nt->Write();
608  ntpad->Write();
609  nthit->Write();
610  ntfinalhit->Write();
611  m_outf->Close();
612  }
614  {
615  EDrift_outf.reset(new TFile("ElectronDriftQA.root", "recreate"));
616  EDrift_outf->cd();
617  deltar->Write();
618  deltaphi->Write();
619  deltaz->Write();
620  deltarnodist->Write();
621  deltaphinodist->Write();
622  deltarnodiff->Write();
623  deltaphinodiff->Write();
624  deltaRphinodiff->Write();
625  deltaphivsRnodiff->Write();
626  hitmapstart->Write();
627  hitmapend->Write();
628  z_startmap->Write();
629  EDrift_outf->Close();
630  }
632 }
633 
634 void PHG4TpcElectronDrift::set_seed(const unsigned int seed)
635 {
636  gsl_rng_set(RandomGenerator.get(), seed);
637 }
638 
640 {
641  // Data on gasses @20 C and 760 Torr from the following source:
642  // http://www.slac.stanford.edu/pubs/icfa/summer98/paper3/paper3.pdf
643  // diffusion and drift velocity for 400kV for NeCF4 50/50 from calculations:
644  // http://skipper.physics.sunysb.edu/~prakhar/tpc/HTML_Gases/split.html
645  static constexpr double Ne_dEdx = 1.56; // keV/cm
646  static constexpr double CF4_dEdx = 7.00; // keV/cm
647  // double Ne_NPrimary = 12; // Number/cm
648  // double CF4_NPrimary = 51; // Number/cm
649  static constexpr double Ne_NTotal = 43; // Number/cm
650  static constexpr double CF4_NTotal = 100; // Number/cm
651  static constexpr double Tpc_NTot = 0.5 * Ne_NTotal + 0.5 * CF4_NTotal;
652  static constexpr double Tpc_dEdx = 0.5 * Ne_dEdx + 0.5 * CF4_dEdx;
653  static constexpr double Tpc_ElectronsPerKeV = Tpc_NTot / Tpc_dEdx;
654  set_default_double_param("diffusion_long", 0.012); // cm/SQRT(cm)
655  set_default_double_param("diffusion_trans", 0.004); // cm/SQRT(cm)
656  set_default_double_param("electrons_per_gev", Tpc_ElectronsPerKeV * 1000000.);
657  set_default_double_param("min_active_radius", 30.); // cm
658  set_default_double_param("max_active_radius", 78.); // cm
659  set_default_double_param("drift_velocity", 8.0 / 1000.0); // cm/ns
660 
661  // These are purely fudge factors, used to increase the resolution to 150 microns and 500 microns, respectively
662  // override them from the macro to get a different resolution
663  set_default_double_param("added_smear_trans", 0.085); // cm
664  set_default_double_param("added_smear_long", 0.105); // cm
665 
666  return;
667 }
668 
670 {
671  m_distortionMap.reset(distortionMap);
672 }
673 
675 {
676  std::cout << "Registering padplane " << std::endl;
677  padplane.reset(inpadplane);
678  padplane->Detector(Detector());
679  padplane->UpdateInternalParameters();
680  std::cout << "padplane registered and parameters updated" << std::endl;
681 
682  return;
683 }