23 #include <phparameter/PHParameterInterface.h>
24 #include <phparameter/PHParameters.h>
25 #include <phparameter/PHParametersContainer.h>
27 #include <pdbcalbase/PdbParameterMapContainer.h>
49 #include <gsl/gsl_randist.h>
50 #include <gsl/gsl_rng.h>
96 std::cout <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
101 const std::string paramnodename =
"G4CELLPARAM_" +
detector;
102 const std::string geonodename =
"G4CELLPAR_" +
detector;
103 const std::string tpcgeonodename =
"G4GEO_" +
detector;
108 std::cout <<
Name() <<
" Could not locate G4HIT node " <<
hitnodename << std::endl;
115 hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
128 DetNode->addNode(newNode);
131 hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode,
"TRKR_HITTRUTHASSOC");
144 DetNode->addNode(newNode);
154 runNode->addNode(newNode);
163 runNode->addNode(RunDetNode);
173 parNode->addNode(ParDetNode);
179 auto tpcparams = findNode::getClass<PHParametersContainer>(ParDetNode, tpcgeonodename);
182 const std::string runparamname =
"G4GEOPARAM_" +
detector;
183 auto tpcpdbparams = findNode::getClass<PdbParameterMapContainer>(RunDetNode, runparamname);
188 tpcparams->CreateAndFillFrom(tpcpdbparams,
detector);
193 std::cout <<
"PHG4TpcElectronDrift::InitRun - failed to find " << runparamname <<
" in order to initialize " << tpcgeonodename <<
". Aborting run ..." << std::endl;
200 const PHParameters *tpcparam = tpcparams->GetParameters(0);
217 se->registerHisto(
dlong);
219 se->registerHisto(
dtrans);
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);
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);
251 padplane->CreateReadoutGeometry(topNode, seggeo);
258 static constexpr
unsigned int print_layer = 18;
267 auto g4hit = findNode::getClass<PHG4HitContainer>(topNode,
hitnodename.c_str());
270 std::cout <<
"Could not locate g4 hit node " <<
hitnodename << std::endl;
276 unsigned int count_g4hits = 0;
277 int count_electrons = 0;
279 double ecollectedhits = 0.0;
280 int ncollectedhits = 0;
282 unsigned int dump_interval = 5000;
283 unsigned int dump_counter = 0;
284 for (
auto hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
289 const double t0 = fmax(hiter->second->get_t(0), hiter->second->get_t(1));
297 double eion = hiter->second->get_eion();
299 count_electrons += n_electrons;
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
313 if (n_electrons == 0)
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;
329 for (
unsigned int i = 0; i < n_electrons; i++)
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));
341 const double rantrans =
347 const double rantime =
350 const double t_final = t_start + t_path + rantime;
351 if (t_final < min_time || t_final >
max_time)
continue;
359 const double radstart = std::sqrt(
square(x_start) +
square(y_start));
360 const double phistart = std::atan2(y_start, x_start);
363 double x_final = x_start + rantrans * std::cos(ranphi);
364 double y_final = y_start + rantrans * std::sin(ranphi);
366 double rad_final = sqrt(
square(x_final) +
square(y_final));
367 double phi_final = atan2(y_final, x_final);
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);
382 rad_final+=r_distortion;
383 phi_final+=phi_distortion;
384 z_final += z_distortion;
386 x_final=rad_final*std::cos(phi_final);
387 y_final=rad_final*std::sin(phi_final);
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);
396 deltaRphinodiff->Fill(radstart, rad_final_nodiff * phi_final_nodiff - radstart * phistart);
401 deltar->Fill(radstart, rad_final - radstart);
402 deltaphi->Fill(phistart, phi_final - phistart);
403 deltaz->Fill(z_start, z_distortion);
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
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;
433 nt->Fill(ihit, t_start, t_final, t_sigma, rad_final, z_start, z_final);
444 single_hitset_iter != single_hitset_range.second;
445 ++single_hitset_iter)
454 std::cout <<
" hitsetkey " << node_hitsetkey <<
" layer " << layer <<
" sector " << sector <<
" side " << side << std::endl;
458 single_hit_iter != single_hit_range.second;
467 std::cout <<
" adding assoc for node_hitsetkey " << node_hitsetkey <<
" single_hitkey " << single_hitkey <<
" g4hitkey " << hiter->first << std::endl;
474 if (dump_counter >= dump_interval || count_g4hits == g4hit->size())
481 temp_hitset_iter != temp_hitset_range.second;
490 std::cout <<
"PHG4TpcElectronDrift: temp_hitset with key: " << node_hitsetkey <<
" in layer " << layer
491 <<
" with sector " << sector <<
" side " << side << std::endl;
499 temp_hit_iter != temp_hit_range.second;
503 TrkrHit *temp_tpchit = temp_hit_iter->second;
504 if (
Verbosity() > 10 && layer == print_layer)
506 std::cout <<
" temp_hitkey " << temp_hitkey <<
" l;ayer " << layer <<
" pad " <<
TpcDefs::getPad(temp_hitkey)
508 <<
" energy " << temp_tpchit->
getEnergy() <<
" eg4hit " << eg4hit << std::endl;
511 ecollectedhits += temp_tpchit->
getEnergy();
516 TrkrHit *node_hit = node_hitsetit->second->getHit(temp_hitkey);
521 node_hitsetit->second->addHitSpecificKey(temp_hitkey, node_hit);
529 if (
Verbosity() > 100 && layer == print_layer)
530 std::cout <<
" ihit " << ihit <<
" collected energy = " << eg4hit << std::endl;
550 std::cout <<
"From PHG4TpcElectronDrift: hitsetcontainer printout at end:" << std::endl;
554 hitset_iter != hitset_range.second;
560 if (layer != print_layer)
continue;
564 std::cout <<
"PHG4TpcElectronDrift: hitset with key: " << hitsetkey <<
" in layer " << layer <<
" with sector " << sector <<
" side " << side << std::endl;
570 hit_iter != hit_range.second;
574 TrkrHit *tpchit = hit_iter->second;
576 <<
" energy " << tpchit->
getEnergy() << std::endl;
583 std::cout <<
"From PHG4TpcElectronDrift: hittruthassoc dump:" << std::endl;
615 EDrift_outf.reset(
new TFile(
"ElectronDriftQA.root",
"recreate"));
645 static constexpr
double Ne_dEdx = 1.56;
646 static constexpr
double CF4_dEdx = 7.00;
649 static constexpr
double Ne_NTotal = 43;
650 static constexpr
double CF4_NTotal = 100;
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;
676 std::cout <<
"Registering padplane " << std::endl;
679 padplane->UpdateInternalParameters();
680 std::cout <<
"padplane registered and parameters updated" << std::endl;