47 #include <gsl/gsl_randist.h>
48 #include <gsl/gsl_rng.h>
51 #include <TMatrixFfwd.h>
88 for(
unsigned int layer = _nlayers_maps + _nlayers_intt;
layer < _nlayers_maps + _nlayers_intt + 16; ++
layer)
98 for(
unsigned int layer = _nlayers_maps + _nlayers_intt +_nlayers_tpc;
layer < _nlayers_maps + _nlayers_intt +_nlayers_tpc + 1; ++
layer)
104 for(
unsigned int layer = _nlayers_maps + _nlayers_intt +_nlayers_tpc + 1;
layer < _nlayers_maps + _nlayers_intt +_nlayers_tpc + 2; ++
layer)
112 for(
unsigned int layer = 0;
layer < _nlayers_maps + _nlayers_intt +_nlayers_tpc + 2; ++
layer)
122 cout <<
"Filling truth cluster node " << endl;
125 auto m_clusterlist = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER_TRUTH");
128 cout <<
PHWHERE <<
" ERROR: Can't find TRKR_CLUSTER_TRUTH" << endl;
135 iter != range.second;
148 bool is_primary =
false;
162 cout <<
PHWHERE <<
" PHG4Particle ID " << gtrackID <<
" gflavor " << gflavor <<
" is_primary " << is_primary
163 <<
" gembed " << gembed << endl;
169 if(
Verbosity() > 0) std::cout <<
" Truth cluster summary for g4particle " << g4particle->
get_track_id() <<
" by layer: " << std::endl;
170 for (
const auto& [ckey, gclus]:truth_clusters )
172 float gx = gclus->getX();
173 float gy = gclus->getY();
174 float gz = gclus->getZ();
175 float ng4hits = gclus->getAdc();
177 TVector3
gpos(gx, gy, gz);
178 float gr = sqrt(gx*gx+gy*gy);
179 float gphi = gpos.Phi();
180 float geta = gpos.Eta();
182 float gphisize = gclus->getSize(1,1);
183 float gzsize = gclus->getSize(2,2);
190 std::cout <<
PHWHERE <<
" **** truth: layer " << layer <<
" truth cluster key " << ckey <<
" ng4hits " << ng4hits << std::endl;
191 std::cout <<
" gr " << gr <<
" gx " << gx <<
" gy " << gy <<
" gz " << gz
192 <<
" gphi " << gphi <<
" geta " << geta <<
" gphisize " << gphisize <<
" gzsize " << gzsize << endl;
198 m_clusterlist->addClusterSpecifyKey(ckey, gclus);
212 for(
auto clusIter = range.first; clusIter != range.second; ++clusIter ){
221 std::cout <<
PHWHERE <<
" copying cluster in layer " << layer <<
" from reco clusters to truth clusters " << std::endl;;
225 m_clusterlist->addClusterSpecifyKey(cluskey, cluster);
231 std::cout <<
"Final TRKR_CLUSTER_TRUTH clusters:";
232 m_clusterlist->identify();
245 std::cout <<
PHWHERE <<
" Truth clustering for particle " << particle->
get_track_id() <<
" with ng4hits " << g4hits.size() << std::endl;;
249 std::map<TrkrDefs::cluskey, TrkrCluster*> truth_clusters;
250 if(g4hits.empty())
return truth_clusters;
263 std::vector<PHG4Hit*> contributing_hits;
264 std::vector<double> contributing_hits_energy;
265 std::vector<std::vector<double>> contributing_hits_entry;
266 std::vector<std::vector<double>> contributing_hits_exit;
267 LayerClusterG4Hits(g4hits, contributing_hits, contributing_hits_energy, contributing_hits_entry, contributing_hits_exit, layer, gx, gy, gz, gt, gedep);
268 if(!(gedep > 0))
continue;
275 unsigned int side = 0;
283 unsigned int stave = 0;
284 unsigned int chip = 0;
285 unsigned int strobe = 0;
291 unsigned int ladderzid = 0;
292 unsigned int ladderphiid = 0;
298 unsigned int tile = 0;
306 std::cout <<
PHWHERE <<
"Bad layer number: " << layer << std::endl;
316 clus->setAdc(adc_output);
317 clus->setPosition(0, gx);
318 clus->setPosition(1, gy);
319 clus->setPosition(2, gz);
331 float g4phisize = NAN;
333 G4ClusterSize(ckey, layer, contributing_hits_entry, contributing_hits_exit, g4phisize, g4zsize);
344 double clusphi = atan2(gy, gx);
351 DIM[1][1] = pow(0.5 * g4phisize, 2);
355 DIM[2][2] = pow(0.5 * g4zsize,2);
369 ROT[0][0] = cos(clusphi);
370 ROT[0][1] = -sin(clusphi);
372 ROT[1][0] = sin(clusphi);
373 ROT[1][1] = cos(clusphi);
379 TMatrixF ROT_T(3, 3);
380 ROT_T.Transpose(ROT);
382 TMatrixF COVAR_DIM(3, 3);
383 COVAR_DIM = ROT * DIM * ROT_T;
385 clus->setSize(0, 0, COVAR_DIM[0][0]);
386 clus->setSize(0, 1, COVAR_DIM[0][1]);
387 clus->setSize(0, 2, COVAR_DIM[0][2]);
388 clus->setSize(1, 0, COVAR_DIM[1][0]);
389 clus->setSize(1, 1, COVAR_DIM[1][1]);
390 clus->setSize(1, 2, COVAR_DIM[1][2]);
391 clus->setSize(2, 0, COVAR_DIM[2][0]);
392 clus->setSize(2, 1, COVAR_DIM[2][1]);
393 clus->setSize(2, 2, COVAR_DIM[2][2]);
396 TMatrixF COVAR_ERR(3, 3);
397 COVAR_ERR = ROT * ERR * ROT_T;
399 clus->setError(0, 0, COVAR_ERR[0][0]);
400 clus->setError(0, 1, COVAR_ERR[0][1]);
401 clus->setError(0, 2, COVAR_ERR[0][2]);
402 clus->setError(1, 0, COVAR_ERR[1][0]);
403 clus->setError(1, 1, COVAR_ERR[1][1]);
404 clus->setError(1, 2, COVAR_ERR[1][2]);
405 clus->setError(2, 0, COVAR_ERR[2][0]);
406 clus->setError(2, 1, COVAR_ERR[2][1]);
407 clus->setError(2, 2, COVAR_ERR[2][2]);
411 std::cout <<
" layer " << layer <<
" cluskey " << ckey <<
" cluster phi " << clusphi <<
" local cluster error rphi " <<
clus_err_rphi[
layer]
415 std::cout <<
" global covariance matrix:" << std::endl;
416 for(
int i1=0;i1<3;++i1)
417 for(
int i2=0;i2<3;++i2)
418 std::cout <<
" " << i1 <<
" " << i2 <<
" cov " << clus->getError(i1,i2) << std::endl;
422 truth_clusters.insert(std::make_pair(ckey, clus));
426 return truth_clusters;
429 void PHTruthClustering::LayerClusterG4Hits(std::set<PHG4Hit*> truth_hits, std::vector<PHG4Hit*> &contributing_hits, std::vector<double> &contributing_hits_energy, std::vector<std::vector<double>> &contributing_hits_entry, std::vector<std::vector<double>> &contributing_hits_exit,
float layer,
float &
x,
float &
y,
float &
z,
float &
t,
float &
e)
457 cout <<
" PHTruthClustering::LayerCluster hits for layer " << layer <<
" with rbin " << rbin <<
" rbout " << rbout << endl;
460 for (std::set<PHG4Hit*>::iterator iter = truth_hits.begin();
461 iter != truth_hits.end();
466 float rbegin = sqrt(this_g4hit->
get_x(0) * this_g4hit->
get_x(0) + this_g4hit->
get_y(0) * this_g4hit->
get_y(0));
467 float rend = sqrt(this_g4hit->
get_x(1) * this_g4hit->
get_x(1) + this_g4hit->
get_y(1) * this_g4hit->
get_y(1));
477 xl[0] = this_g4hit->
get_x(0);
478 yl[0] = this_g4hit->
get_y(0);
479 zl[0] = this_g4hit->
get_z(0);
480 xl[1] = this_g4hit->
get_x(1);
481 yl[1] = this_g4hit->
get_y(1);
482 zl[1] = this_g4hit->
get_z(1);
486 xl[0] = this_g4hit->
get_x(1);
487 yl[0] = this_g4hit->
get_y(1);
488 zl[0] = this_g4hit->
get_z(1);
489 xl[1] = this_g4hit->
get_x(0);
490 yl[1] = this_g4hit->
get_y(0);
491 zl[1] = this_g4hit->
get_z(0);
497 if (rbegin < rbin && rend < rbin)
499 if (rbegin > rbout && rend > rbout)
505 cout <<
" keep g4hit with rbegin " << rbegin <<
" rend " << rend
506 <<
" xbegin " << xl[0] <<
" xend " << xl[1]
507 <<
" ybegin " << yl[0] <<
" yend " << yl[1]
508 <<
" zbegin " << zl[0] <<
" zend " << zl[1]
528 xin = xl[0] + t * (xl[1] - xl[0]);
529 yin = yl[0] + t * (yl[1] - yl[0]);
530 zin = zl[0] + t * (zl[1] - zl[0]);
540 xout = xl[0] + t * (xl[1] - xl[0]);
541 yout = yl[0] + t * (yl[1] - yl[0]);
542 zout = zl[0] + t * (zl[1] - zl[0]);
546 double rin = sqrt(xin*xin + yin*yin);
547 double rout = sqrt(xout*xout + yout*yout);
550 double efrac = this_g4hit->
get_edep() * (rout - rin) / (rend - rbegin);
551 gx += (xin + xout) * 0.5 * efrac;
552 gy += (yin + yout) * 0.5 * efrac;
553 gz += (zin + zout) * 0.5 * efrac;
555 gr += (rin + rout) * 0.5 * efrac;
560 cout <<
" rin " << rin <<
" rout " << rout
561 <<
" xin " << xin <<
" xout " << xout <<
" yin " << yin <<
" yout " << yout <<
" zin " << zin <<
" zout " << zout
562 <<
" edep " << this_g4hit->
get_edep()
563 <<
" this_edep " << efrac << endl;
564 cout <<
" xavge " << (xin+xout) * 0.5 <<
" yavge " << (yin+yout) * 0.5 <<
" zavge " << (zin+zout) * 0.5 <<
" ravge " << (rin+rout) * 0.5
569 std::vector<double> entry_loc;
570 entry_loc.push_back(xin);
571 entry_loc.push_back(yin);
572 entry_loc.push_back(zin);
573 std::vector<double> exit_loc;
574 exit_loc.push_back(xout);
575 exit_loc.push_back(yout);
576 exit_loc.push_back(zout);
579 contributing_hits.push_back(this_g4hit);
580 contributing_hits_energy.push_back( this_g4hit->
get_edep() * (zout - zin) / (zl[1] - zl[0]) );
581 contributing_hits_entry.push_back(entry_loc);
582 contributing_hits_exit.push_back(exit_loc);
595 gr = (rbin + rbout) * 0.5;
600 cout <<
" weighted means: gx " << gx <<
" gy " << gy <<
" gz " << gz <<
" gr " << gr <<
" e " << gwt << endl;
607 float rentry = 999.0;
608 float xentry = 999.0;
609 float yentry = 999.0;
610 float zentry = 999.0;
611 float rexit = - 999.0;
612 float xexit = -999.0;
613 float yexit = -999.0;
614 float zexit = -999.0;
616 for(
unsigned int ientry = 0; ientry < contributing_hits_entry.size(); ++ientry)
618 float tmpx = contributing_hits_entry[ientry][0];
619 float tmpy = contributing_hits_entry[ientry][1];
620 float tmpr = sqrt(tmpx*tmpx + tmpy*tmpy);
625 xentry = contributing_hits_entry[ientry][0];
626 yentry = contributing_hits_entry[ientry][1];
627 zentry = contributing_hits_entry[ientry][2];
630 tmpx = contributing_hits_exit[ientry][0];
631 tmpy = contributing_hits_exit[ientry][1];
632 tmpr = sqrt(tmpx*tmpx + tmpy*tmpy);
637 xexit = contributing_hits_exit[ientry][0];
638 yexit = contributing_hits_exit[ientry][1];
639 zexit = contributing_hits_exit[ientry][2];
643 float geo_r = (rentry+rexit)*0.5;
644 float geo_x = (xentry+xexit)*0.5;
645 float geo_y = (yentry+yexit)*0.5;
646 float geo_z = (zentry+zexit)*0.5;
659 cout <<
" rentry " << rentry <<
" rexit " << rexit
660 <<
" xentry " << xentry <<
" xexit " << xexit <<
" yentry " << yentry <<
" yexit " << yexit <<
" zentry " << zentry <<
" zexit " << zexit << endl;
662 cout <<
" geometric means: geo_x " << geo_x <<
" geo_y " << geo_y <<
" geo_z " << geo_z <<
" geo r " << geo_r <<
" e " << gwt << endl << endl;
670 for (std::set<PHG4Hit*>::iterator iter = truth_hits.begin();
671 iter != truth_hits.end();
677 if(this_g4hit->
get_layer() != (
unsigned int) layer)
continue;
686 std::vector<double> entry_loc;
687 entry_loc.push_back(this_g4hit->
get_x(0));
688 entry_loc.push_back(this_g4hit->
get_y(0));
689 entry_loc.push_back(this_g4hit->
get_z(0));
690 std::vector<double> exit_loc;
691 exit_loc.push_back(this_g4hit->
get_x(1));
692 exit_loc.push_back(this_g4hit->
get_y(1));
693 exit_loc.push_back(this_g4hit->
get_z(1));
696 contributing_hits.push_back(this_g4hit);
697 contributing_hits_energy.push_back( this_g4hit->
get_edep() );
698 contributing_hits_entry.push_back(entry_loc);
699 contributing_hits_exit.push_back(exit_loc);
717 double inner_x = NAN;
718 double inner_y = NAN;
719 double inner_z = NAN;;
722 double outer_x = NAN;
723 double outer_y = NAN;
724 double outer_z = NAN;
726 for(
unsigned int ihit=0;
ihit<contributing_hits_entry.size(); ++
ihit)
728 double rad1 = sqrt(pow(contributing_hits_entry[
ihit][0], 2) + pow(contributing_hits_entry[ihit][1], 2));
729 if(rad1 < inner_radius)
732 inner_x = contributing_hits_entry[
ihit][0];
733 inner_y = contributing_hits_entry[
ihit][1];
734 inner_z = contributing_hits_entry[
ihit][2];
737 double rad2 = sqrt(pow(contributing_hits_exit[ihit][0], 2) + pow(contributing_hits_exit[ihit][1], 2));
738 if(rad2 > outer_radius)
741 outer_x = contributing_hits_exit[
ihit][0];
742 outer_y = contributing_hits_exit[
ihit][1];
743 outer_z = contributing_hits_exit[
ihit][2];
747 double inner_phi = atan2(inner_y, inner_x);
748 double outer_phi = atan2(outer_y, outer_x);
749 double avge_z = (outer_z + inner_z) / 2.0;
756 if(radius > 28 && radius < 80)
760 double tpc_length = 211.0;
761 double drift_velocity = 8.0 / 1000.0;
765 double diffusion_trans = 0.006;
766 double phidiffusion = diffusion_trans * sqrt(tpc_length / 2. - fabs(avge_z));
768 double added_smear_trans = 0.085;
769 double gem_spread = 0.04;
771 if(outer_phi < inner_phi)
swap(outer_phi, inner_phi);
774 double g4max_phi = outer_phi + sigmas * sqrt( pow(phidiffusion, 2) + pow(added_smear_trans, 2) + pow(gem_spread, 2) ) / radius;
775 double g4min_phi = inner_phi - sigmas * sqrt( pow(phidiffusion, 2) + pow(added_smear_trans, 2) + pow(gem_spread, 2) ) / radius;
778 unsigned int phibinmin = layergeom->
get_phibin(g4min_phi);
779 unsigned int phibinmax = layergeom->
get_phibin(g4max_phi);
780 unsigned int phibinwidth = phibinmax - phibinmin + 1;
788 outer_z = fabs(outer_z);
789 inner_z = fabs(inner_z);
791 double diffusion_long = 0.015;
792 double zdiffusion = diffusion_long * sqrt(tpc_length / 2. - fabs(avge_z)) ;
793 double zshaping_lead = 32.0 * drift_velocity;
794 double zshaping_tail = 48.0 * drift_velocity;
795 double added_smear_long = 0.105;
798 if(outer_z < inner_z)
swap(outer_z, inner_z);
799 g4max_z = outer_z + sigmas*sqrt(pow(zdiffusion,2) + pow(added_smear_long,2) + pow(zshaping_lead, 2));
800 g4min_z = inner_z - sigmas*sqrt(pow(zdiffusion,2) + pow(added_smear_long,2) + pow(zshaping_tail, 2));
803 unsigned int binmin = layergeom->
get_zbin(g4min_z);
804 unsigned int binmax = layergeom->
get_zbin(g4max_z);
805 if(binmax < binmin)
swap(binmax, binmin);
806 unsigned int binwidth = binmax - binmin + 1;
811 else if(radius > 5 && radius < 20)
818 double world_inner[3] = {inner_x, inner_y, inner_z};
819 TVector3 world_inner_vec = {inner_x, inner_y, inner_z};
821 int segment_z_bin, segment_phi_bin;
826 double yin = local_inner_vec[1];
827 double zin = local_inner_vec[2];
828 int strip_y_index, strip_z_index;
832 double world_outer[3] = {outer_x, outer_y, outer_z};
833 TVector3 world_outer_vec = {outer_x, outer_y, outer_z};
839 double yout = local_outer_vec[1];
840 double zout = local_outer_vec[2];
841 int strip_y_index_out, strip_z_index_out;
844 int strips = abs(strip_y_index_out - strip_y_index) + 1;
845 int cols = abs(strip_z_index_out - strip_z_index) + 1;
851 g4phisize = strip_width;
852 g4zsize = strip_length;
868 unsigned int stave, stave_outer;
869 unsigned int chip, chip_outer;
871 int column, column_outer;
874 double max_diffusion_radius = 25.0e-4;
875 double min_diffusion_radius = 8.0e-4;
879 TVector3 world_inner = {inner_x, inner_y, inner_z};
880 std::vector<double> world_inner_vec = { world_inner[0], world_inner[1], world_inner[2] };
886 TVector3 world_outer = {outer_x, outer_y, outer_z};
887 std::vector<double> world_outer_vec = { world_outer[0], world_outer[1], world_outer[2] };
893 double diff = max_diffusion_radius * 0.6;
894 if(local_outer[0] < local_inner[0])
896 local_outer[0] += diff;
897 local_inner[0] -= diff;
899 double diff_outer = min_diffusion_radius * 0.6;
900 if(local_outer[2] < local_inner[2])
901 diff_outer = -diff_outer;
902 local_outer[2] += diff_outer;
903 local_inner[2] -= diff_outer;
908 if(row_outer < row)
swap(row_outer, row);
909 unsigned int rows = row_outer - row + 1;
912 if(column_outer < column)
swap(column_outer, column);
913 unsigned int columns = column_outer - column + 1;
927 std::set<PHG4Hit*> truth_hits;
936 PHG4Hit* g4hit = g4iter->second;
939 truth_hits.insert(g4hit);
951 PHG4Hit* g4hit = g4iter->second;
954 truth_hits.insert(g4hit);
966 PHG4Hit* g4hit = g4iter->second;
969 truth_hits.insert(g4hit);
981 PHG4Hit* g4hit = g4iter->second;
984 truth_hits.insert(g4hit);
1002 float A = (x[1] - x[0]) * (x[1] - x[0]) + (y[1] - y[0]) * (y[1] - y[0]);
1003 float B = 2.0 * x[0] * (x[1] - x[0]) + 2.0 * y[0] * (y[1] - y[0]);
1004 float C = x[0] * x[0] + y[0] * y[0] - radius * radius;
1005 float tup = (-B + sqrt(B * B - 4.0 * A * C)) / (2.0 * A);
1006 float tdn = (-B - sqrt(B * B - 4.0 * A * C)) / (2.0 * A);
1010 if (tdn >= -0.0
e-4 && tdn <= 1.0004)
1012 else if (tup >= -0.0
e-4 && tup <= 1.0004)
1016 cout <<
PHWHERE <<
" **** Oops! No valid solution for tup or tdn, tdn = " << tdn <<
" tup = " << tup << endl;
1017 cout <<
" radius " << radius <<
" rbegin " << sqrt(x[0] * x[0] + y[0] * y[0]) <<
" rend " << sqrt(x[1] * x[1] + y[1] * y[1]) << endl;
1018 cout <<
" x0 " << x[0] <<
" x1 " << x[1] << endl;
1019 cout <<
" y0 " << y[0] <<
" y1 " << y[1] << endl;
1020 cout <<
" z0 " << z[0] <<
" z1 " << z[1] << endl;
1021 cout <<
" A " << A <<
" B " << B <<
" C " << C << endl;
1030 double phi = atan2(y, x);
1031 unsigned int sector = (int) (12.0 * (phi + M_PI) / (2.0 * M_PI) );
1044 double Tpc_NTot = 0.5*Ne_NTotal + 0.5*
CF4_NTotal;
1045 double Tpc_dEdx = 0.5*Ne_dEdx + 0.5*
CF4_dEdx;
1046 double Tpc_ElectronsPerKeV = Tpc_NTot / Tpc_dEdx;
1047 double electrons_per_gev = Tpc_ElectronsPerKeV * 1e6;
1049 double gem_amplification = 1400;
1050 double input_electrons = gedep * electrons_per_gev * gem_amplification;
1053 double ChargeToPeakVolts = 20;
1054 double ADCSignalConversionGain = ChargeToPeakVolts * 1.60e-04 * 2.4;
1055 double adc_input_voltage = input_electrons * ADCSignalConversionGain;
1056 unsigned int adc_output = (
unsigned int) (adc_input_voltage * 1024.0 / 2200.0);
1057 if (adc_output > 1023) adc_output = 1023;
1064 _tgeometry = findNode::getClass<ActsGeometry>(topNode,
"ActsGeometry");
1065 _g4truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode,
"G4TruthInfo");
1066 if (!_g4truth_container)
1068 cerr <<
PHWHERE <<
" ERROR: Can't find node G4TruthInfo" << endl;
1072 _g4hits_mms = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_MICROMEGAS");
1073 _g4hits_svtx = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_TPC");
1074 _g4hits_tracker = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_INTT");
1075 _g4hits_maps = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_MVTX");
1077 _mms_geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode,
"CYLINDERGEOM_MICROMEGAS_FULL");
1078 _tpc_geom_container = findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode,
"CYLINDERCELLGEOM_SVTX");
1079 _intt_geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode,
"CYLINDERGEOM_INTT");
1080 _mvtx_geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode,
"CYLINDERGEOM_MVTX");
1088 cout <<
PHWHERE <<
"DST Node missing, doing nothing." << endl;
1092 auto trkrclusters = findNode::getClass<TrkrClusterContainer>(dstNode,
"TRKR_CLUSTER_TRUTH");
1101 dstNode->addNode(DetNode);
1107 DetNode->
addNode(TrkrClusterContainerNode);
1110 _reco_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
1113 cerr <<
PHWHERE <<
" ERROR: Can't find node TRKR_CLUSTER" << endl;