45 : _basetrutheval(topNode)
51 , _cache_all_truth_hits()
52 , _cache_all_truth_hits_g4particle()
53 , _cache_all_truth_clusters_g4particle()
54 , _cache_get_innermost_truth_hit()
55 , _cache_get_outermost_truth_hit()
56 , _cache_get_primary_particle_g4hit()
68 cout <<
"SvtxTruthEval::~SvtxTruthEval() - Error Count: " <<
_errors << endl;
93 return std::set<PHG4Hit*>();
107 std::set<PHG4Hit*> truth_hits;
116 PHG4Hit* g4hit = g4iter->second;
117 truth_hits.insert(g4hit);
128 PHG4Hit* g4hit = g4iter->second;
129 truth_hits.insert(g4hit);
140 PHG4Hit* g4hit = g4iter->second;
141 truth_hits.insert(g4hit);
152 PHG4Hit* g4hit = g4iter->second;
153 truth_hits.insert(g4hit);
167 return std::set<PHG4Hit*>();
177 return std::set<PHG4Hit*>();
186 std::map<PHG4Particle*, std::set<PHG4Hit*> >::iterator iter =
195 std::set<PHG4Hit*> truth_hits;
201 std::multimap<const int, PHG4Hit*> temp_clusters_from_particles;
208 PHG4Hit* g4hit = g4iter->second;
209 temp_clusters_from_particles.insert(make_pair(g4hit->
get_trkid(),g4hit));
220 PHG4Hit* g4hit = g4iter->second;
221 temp_clusters_from_particles.insert(make_pair(g4hit->
get_trkid(),g4hit));
232 PHG4Hit* g4hit = g4iter->second;
233 temp_clusters_from_particles.insert(make_pair(g4hit->
get_trkid(),g4hit));
244 PHG4Hit* g4hit = g4iter->second;
245 temp_clusters_from_particles.insert(make_pair(g4hit->
get_trkid(),g4hit));
252 iter != range.second; ++iter){
254 std::set<PHG4Hit*> truth_hits;
255 std::multimap<const int, PHG4Hit*>::const_iterator lower_bound = temp_clusters_from_particles.lower_bound(g4particle->
get_track_id());
256 std::multimap<const int, PHG4Hit*>::const_iterator upper_bound = temp_clusters_from_particles.upper_bound(g4particle->
get_track_id());
257 std::multimap<const int, PHG4Hit*>::const_iterator cfp_iter;
258 for(cfp_iter = lower_bound;cfp_iter != upper_bound;++cfp_iter){
259 PHG4Hit* g4hit = cfp_iter->second;
260 truth_hits.insert(g4hit);
272 return std::map<unsigned int, std::shared_ptr<TrkrCluster> >();
282 return std::map<unsigned int, std::shared_ptr<TrkrCluster> >();
287 std::map<PHG4Particle*, std::map<unsigned int, std::shared_ptr<TrkrCluster> > >::iterator iter =
301 float ng4hits = g4hits.size();
303 return std::map<unsigned int, std::shared_ptr<TrkrCluster> >();
307 std::map<unsigned int, std::shared_ptr<TrkrCluster>> truth_clusters;
321 std::vector<PHG4Hit*> contributing_hits;
322 std::vector<double> contributing_hits_energy;
323 std::vector<std::vector<double>> contributing_hits_entry;
324 std::vector<std::vector<double>> contributing_hits_exit;
325 LayerClusterG4Hits(g4hits, contributing_hits, contributing_hits_energy, contributing_hits_entry, contributing_hits_exit, layer, gx, gy, gz, gt, gedep);
326 if(!(gedep > 0))
continue;
333 unsigned int side = 0;
336 unsigned int sector = 0;
341 unsigned int stave = 0;
342 unsigned int chip = 0;
348 unsigned int ladderzid = 0;
349 unsigned int ladderphiid = 0;
354 unsigned int tile = 0;
362 std::cout <<
PHWHERE <<
"Bad layer number: " << layer << std::endl;
367 clus->setClusKey(ckey);
373 clus->setAdc(adc_value);
374 clus->setPosition(0, gx);
375 clus->setPosition(1, gy);
376 clus->setPosition(2, gz);
380 for(
unsigned int i=0; i< contributing_hits.size(); ++i)
386 float g4phisize = NAN;
388 G4ClusterSize(layer, contributing_hits_entry, contributing_hits_exit, g4phisize, g4zsize);
390 for(
int i1=0;i1<3;++i1)
391 for(
int i2=0;i2<3;++i2)
393 clus->setSize(i1, i2, 0.0);
394 clus->setError(i1, i2, 0.0);
396 clus->setError(0,0,gedep);
397 clus->setSize(1, 1, g4phisize);
398 clus->setSize(2, 2, g4zsize);
399 clus->setError(1, 1, g4phisize/sqrt(12));
400 clus->setError(2, 2, g4zsize/sqrt(12.0));
402 truth_clusters.insert(make_pair(layer, clus));
408 return truth_clusters;
411 void SvtxTruthEval::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)
439 cout <<
" TruthEval::LayerCluster hits for layer " << layer <<
" with rbin " << rbin <<
" rbout " << rbout << endl;
442 for (std::set<PHG4Hit*>::iterator iter = truth_hits.begin();
443 iter != truth_hits.end();
448 float rbegin = sqrt(this_g4hit->
get_x(0) * this_g4hit->
get_x(0) + this_g4hit->
get_y(0) * this_g4hit->
get_y(0));
449 float rend = sqrt(this_g4hit->
get_x(1) * this_g4hit->
get_x(1) + this_g4hit->
get_y(1) * this_g4hit->
get_y(1));
459 xl[0] = this_g4hit->
get_x(0);
460 yl[0] = this_g4hit->
get_y(0);
461 zl[0] = this_g4hit->
get_z(0);
462 xl[1] = this_g4hit->
get_x(1);
463 yl[1] = this_g4hit->
get_y(1);
464 zl[1] = this_g4hit->
get_z(1);
468 xl[0] = this_g4hit->
get_x(1);
469 yl[0] = this_g4hit->
get_y(1);
470 zl[0] = this_g4hit->
get_z(1);
471 xl[1] = this_g4hit->
get_x(0);
472 yl[1] = this_g4hit->
get_y(0);
473 zl[1] = this_g4hit->
get_z(0);
479 if (rbegin < rbin && rend < rbin)
481 if (rbegin > rbout && rend > rbout)
487 cout <<
" keep g4hit with rbegin " << rbegin <<
" rend " << rend
488 <<
" xbegin " << xl[0] <<
" xend " << xl[1]
489 <<
" ybegin " << yl[0] <<
" yend " << yl[1]
490 <<
" zbegin " << zl[0] <<
" zend " << zl[1]
510 xin = xl[0] + t * (xl[1] - xl[0]);
511 yin = yl[0] + t * (yl[1] - yl[0]);
512 zin = zl[0] + t * (zl[1] - zl[0]);
522 xout = xl[0] + t * (xl[1] - xl[0]);
523 yout = yl[0] + t * (yl[1] - yl[0]);
524 zout = zl[0] + t * (zl[1] - zl[0]);
528 double rin = sqrt(xin*xin + yin*yin);
529 double rout = sqrt(xout*xout + yout*yout);
532 double efrac = this_g4hit->
get_edep() * (rout - rin) / (rend - rbegin);
533 gx += (xin + xout) * 0.5 * efrac;
534 gy += (yin + yout) * 0.5 * efrac;
535 gz += (zin + zout) * 0.5 * efrac;
537 gr += (rin + rout) * 0.5 * efrac;
542 cout <<
" rin " << rin <<
" rout " << rout
543 <<
" xin " << xin <<
" xout " << xout <<
" yin " << yin <<
" yout " << yout <<
" zin " << zin <<
" zout " << zout
544 <<
" edep " << this_g4hit->
get_edep()
545 <<
" this_edep " << efrac << endl;
546 cout <<
" xavge " << (xin+xout) * 0.5 <<
" yavge " << (yin+yout) * 0.5 <<
" zavge " << (zin+zout) * 0.5 <<
" ravge " << (rin+rout) * 0.5
551 std::vector<double> entry_loc;
552 entry_loc.push_back(xin);
553 entry_loc.push_back(yin);
554 entry_loc.push_back(zin);
555 std::vector<double> exit_loc;
556 exit_loc.push_back(xout);
557 exit_loc.push_back(yout);
558 exit_loc.push_back(zout);
561 contributing_hits.push_back(this_g4hit);
562 contributing_hits_energy.push_back( this_g4hit->
get_edep() * (zout - zin) / (zl[1] - zl[0]) );
563 contributing_hits_entry.push_back(entry_loc);
564 contributing_hits_exit.push_back(exit_loc);
577 gr = (rbin + rbout) * 0.5;
582 cout <<
" weighted means: gx " << gx <<
" gy " << gy <<
" gz " << gz <<
" gr " << gr <<
" e " << gwt << endl;
589 float rentry = 999.0;
590 float xentry = 999.0;
591 float yentry = 999.0;
592 float zentry = 999.0;
593 float rexit = - 999.0;
594 float xexit = -999.0;
595 float yexit = -999.0;
596 float zexit = -999.0;
598 for(
unsigned int ientry = 0; ientry < contributing_hits_entry.size(); ++ientry)
600 float tmpx = contributing_hits_entry[ientry][0];
601 float tmpy = contributing_hits_entry[ientry][1];
602 float tmpr = sqrt(tmpx*tmpx + tmpy*tmpy);
607 xentry = contributing_hits_entry[ientry][0];
608 yentry = contributing_hits_entry[ientry][1];
609 zentry = contributing_hits_entry[ientry][2];
612 tmpx = contributing_hits_exit[ientry][0];
613 tmpy = contributing_hits_exit[ientry][1];
614 tmpr = sqrt(tmpx*tmpx + tmpy*tmpy);
619 xexit = contributing_hits_exit[ientry][0];
620 yexit = contributing_hits_exit[ientry][1];
621 zexit = contributing_hits_exit[ientry][2];
625 float geo_r = (rentry+rexit)*0.5;
626 float geo_x = (xentry+xexit)*0.5;
627 float geo_y = (yentry+yexit)*0.5;
628 float geo_z = (zentry+zexit)*0.5;
641 cout <<
" rentry " << rentry <<
" rexit " << rexit
642 <<
" xentry " << xentry <<
" xexit " << xexit <<
" yentry " << yentry <<
" yexit " << yexit <<
" zentry " << zentry <<
" zexit " << zexit << endl;
644 cout <<
" geometric means: geo_x " << geo_x <<
" geo_y " << geo_y <<
" geo_z " << geo_z <<
" geo r " << geo_r <<
" e " << gwt << endl << endl;
652 for (std::set<PHG4Hit*>::iterator iter = truth_hits.begin();
653 iter != truth_hits.end();
659 if(this_g4hit->
get_layer() != (
unsigned int) layer)
continue;
668 std::vector<double> entry_loc;
669 entry_loc.push_back(this_g4hit->
get_x(0));
670 entry_loc.push_back(this_g4hit->
get_y(0));
671 entry_loc.push_back(this_g4hit->
get_z(0));
672 std::vector<double> exit_loc;
673 exit_loc.push_back(this_g4hit->
get_x(1));
674 exit_loc.push_back(this_g4hit->
get_y(1));
675 exit_loc.push_back(this_g4hit->
get_z(1));
678 contributing_hits.push_back(this_g4hit);
679 contributing_hits_energy.push_back( this_g4hit->
get_edep() );
680 contributing_hits_entry.push_back(entry_loc);
681 contributing_hits_exit.push_back(exit_loc);
694 void SvtxTruthEval::G4ClusterSize(
unsigned int layer, std::vector<std::vector<double>> contributing_hits_entry,std::vector<std::vector<double>> contributing_hits_exit,
float &g4phisize,
float &g4zsize)
699 double inner_x = NAN;
700 double inner_y = NAN;
701 double inner_z = NAN;;
704 double outer_x = NAN;
705 double outer_y = NAN;
706 double outer_z = NAN;
708 for(
unsigned int ihit=0;
ihit<contributing_hits_entry.size(); ++
ihit)
710 double rad1 = sqrt(pow(contributing_hits_entry[
ihit][0], 2) + pow(contributing_hits_entry[ihit][1], 2));
711 if(rad1 < inner_radius)
714 inner_x = contributing_hits_entry[
ihit][0];
715 inner_y = contributing_hits_entry[
ihit][1];
716 inner_z = contributing_hits_entry[
ihit][2];
719 double rad2 = sqrt(pow(contributing_hits_exit[ihit][0], 2) + pow(contributing_hits_exit[ihit][1], 2));
720 if(rad2 > outer_radius)
723 outer_x = contributing_hits_exit[
ihit][0];
724 outer_y = contributing_hits_exit[
ihit][1];
725 outer_z = contributing_hits_exit[
ihit][2];
729 double inner_phi = atan2(inner_y, inner_x);
730 double outer_phi = atan2(outer_y, outer_x);
731 double avge_z = (outer_z + inner_z) / 2.0;
738 if(radius > 28 && radius < 80)
742 double tpc_length = 211.0;
743 double drift_velocity = 8.0 / 1000.0;
747 double diffusion_trans = 0.006;
748 double phidiffusion = diffusion_trans * sqrt(tpc_length / 2. - fabs(avge_z));
750 double added_smear_trans = 0.085;
751 double gem_spread = 0.04;
753 if(outer_phi < inner_phi)
swap(outer_phi, inner_phi);
756 double g4max_phi = outer_phi + sigmas * sqrt( pow(phidiffusion, 2) + pow(added_smear_trans, 2) + pow(gem_spread, 2) ) /
radius;
757 double g4min_phi = inner_phi - sigmas * sqrt( pow(phidiffusion, 2) + pow(added_smear_trans, 2) + pow(gem_spread, 2) ) /
radius;
760 unsigned int phibinmin = layergeom->
get_phibin(g4min_phi);
761 unsigned int phibinmax = layergeom->
get_phibin(g4max_phi);
762 unsigned int phibinwidth = phibinmax - phibinmin + 1;
770 outer_z = fabs(outer_z);
771 inner_z = fabs(inner_z);
773 double diffusion_long = 0.015;
774 double zdiffusion = diffusion_long * sqrt(tpc_length / 2. - fabs(avge_z)) ;
775 double zshaping_lead = 32.0 * drift_velocity;
776 double zshaping_tail = 48.0 * drift_velocity;
777 double added_smear_long = 0.105;
780 if(outer_z < inner_z)
swap(outer_z, inner_z);
781 g4max_z = outer_z + sigmas*sqrt(pow(zdiffusion,2) + pow(added_smear_long,2) + pow(zshaping_lead, 2));
782 g4min_z = inner_z - sigmas*sqrt(pow(zdiffusion,2) + pow(added_smear_long,2) + pow(zshaping_tail, 2));
785 unsigned int binmin = layergeom->
get_zbin(g4min_z);
786 unsigned int binmax = layergeom->
get_zbin(g4max_z);
787 if(binmax < binmin)
swap(binmax, binmin);
788 unsigned int binwidth = binmax - binmin + 1;
791 g4zsize = (double) binwidth * layergeom->
get_zstep();
793 else if(radius > 5 && radius < 20)
800 double world_inner[3] = {inner_x, inner_y, inner_z};
801 TVector3 world_inner_vec = {inner_x, inner_y, inner_z};
803 int segment_z_bin, segment_phi_bin;
807 double yin = local_inner_vec[1];
808 double zin = local_inner_vec[2];
809 int strip_y_index, strip_z_index;
813 double world_outer[3] = {outer_x, outer_y, outer_z};
814 TVector3 world_outer_vec = {outer_x, outer_y, outer_z};
819 double yout = local_outer_vec[1];
820 double zout = local_outer_vec[2];
821 int strip_y_index_out, strip_z_index_out;
824 int strips =
abs(strip_y_index_out - strip_y_index) + 1;
825 int cols =
abs(strip_z_index_out - strip_z_index) + 1;
831 g4phisize = strip_width;
832 g4zsize = strip_length;
848 unsigned int stave, stave_outer;
849 unsigned int chip, chip_outer;
851 int column, column_outer;
854 double max_diffusion_radius = 25.0e-4;
855 double min_diffusion_radius = 8.0e-4;
859 TVector3 world_inner = {inner_x, inner_y, inner_z};
860 std::vector<double> world_inner_vec = { world_inner[0], world_inner[1], world_inner[2] };
864 TVector3 world_outer = {outer_x, outer_y, outer_z};
865 std::vector<double> world_outer_vec = { world_outer[0], world_outer[1], world_outer[2] };
869 double diff = max_diffusion_radius * 0.6;
870 if(local_outer[0] < local_inner[0])
872 local_outer[0] += diff;
873 local_inner[0] -= diff;
875 double diff_outer = min_diffusion_radius * 0.6;
876 if(local_outer[2] < local_inner[2])
877 diff_outer = -diff_outer;
878 local_outer[2] += diff_outer;
879 local_inner[2] -= diff_outer;
884 if(row_outer < row)
swap(row_outer, row);
885 unsigned int rows = row_outer - row + 1;
886 g4phisize = (double) rows * layergeom->
get_pixel_x();
888 if(column_outer < column)
swap(column_outer, column);
889 unsigned int columns = column_outer - column + 1;
890 g4zsize = (double) columns * layergeom->
get_pixel_z();
904 std::set<PHG4Hit *> g4hit_set;
907 std::pair<std::multimap<TrkrDefs::cluskey, PHG4Hit*>::iterator,
910 for(std::multimap<TrkrDefs::cluskey, PHG4Hit*>::iterator jter = ret.first; jter != ret.second; ++jter)
912 g4hit_set.insert(jter->second);
936 PHG4Hit* innermost_hit =
nullptr;
937 float innermost_radius =
FLT_MAX;
940 for (std::set<PHG4Hit*>::iterator iter = truth_hits.begin();
941 iter != truth_hits.end();
945 float x = candidate->
get_x(0);
946 float y = candidate->
get_y(0);
947 float r = sqrt(x * x + y * y);
948 if (r < innermost_radius)
950 innermost_radius =
r;
951 innermost_hit = candidate;
955 return innermost_hit;
976 PHG4Hit* outermost_hit =
nullptr;
977 float outermost_radius =
FLT_MAX * -1.0;
981 std::map<PHG4Particle*, PHG4Hit*>::iterator iter =
990 for (std::set<PHG4Hit*>::iterator iter = truth_hits.begin();
991 iter != truth_hits.end();
995 float x = candidate->
get_x(1);
996 float y = candidate->
get_y(1);
997 float r = sqrt(x * x + y * y);
998 if (r > outermost_radius)
1000 outermost_radius =
r;
1001 outermost_hit = candidate;
1006 return outermost_hit;
1049 std::map<PHG4Hit*, PHG4Particle*>::iterator iter =
1053 return iter->second;
1100 _truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,
"G4TruthInfo");
1102 _g4hits_mms = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_MICROMEGAS");
1103 _g4hits_svtx = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_TPC");
1104 _g4hits_tracker = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_INTT");
1105 _g4hits_maps = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_MVTX");
1107 _mms_geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode,
"CYLINDERGEOM_MICROMEGAS_FULL");
1108 _tpc_geom_container = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode,
"CYLINDERCELLGEOM_SVTX");
1109 _intt_geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode,
"CYLINDERGEOM_INTT");
1110 _mvtx_geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode,
"CYLINDERGEOM_MVTX");
1119 else if (!_truthinfo)
1140 float A = (x[1] - x[0]) * (x[1] - x[0]) + (y[1] - y[0]) * (y[1] - y[0]);
1141 float B = 2.0 * x[0] * (x[1] - x[0]) + 2.0 * y[0] * (y[1] - y[0]);
1142 float C = x[0] * x[0] + y[0] * y[0] - radius *
radius;
1143 float tup = (-B + sqrt(B * B - 4.0 * A * C)) / (2.0 * A);
1144 float tdn = (-B - sqrt(B * B - 4.0 * A * C)) / (2.0 * A);
1148 if (tdn >= -0.0
e-4 && tdn <= 1.0004)
1150 else if (tup >= -0.0
e-4 && tup <= 1.0004)
1154 cout <<
PHWHERE <<
" **** Oops! No valid solution for tup or tdn, tdn = " << tdn <<
" tup = " << tup << endl;
1155 cout <<
" radius " << radius <<
" rbegin " << sqrt(x[0] * x[0] + y[0] * y[0]) <<
" rend " << sqrt(x[1] * x[1] + y[1] * y[1]) << endl;
1156 cout <<
" x0 " << x[0] <<
" x1 " << x[1] << endl;
1157 cout <<
" y0 " << y[0] <<
" y1 " << y[1] << endl;
1158 cout <<
" z0 " << z[0] <<
" z1 " << z[1] << endl;
1159 cout <<
" A " << A <<
" B " << B <<
" C " << C << endl;
1172 double Ne_dEdx = 1.56;
1173 double CF4_dEdx = 7.00;
1174 double Ne_NTotal = 43;
1175 double CF4_NTotal = 100;
1176 double Tpc_NTot = 0.5*Ne_NTotal + 0.5*CF4_NTotal;
1177 double Tpc_dEdx = 0.5*Ne_dEdx + 0.5*CF4_dEdx;
1178 double Tpc_ElectronsPerKeV = Tpc_NTot / Tpc_dEdx;
1179 double electrons_per_gev = Tpc_ElectronsPerKeV * 1
e6;
1181 double gem_amplification = 1400;
1182 double input_electrons = gedep * electrons_per_gev * gem_amplification;
1185 double ChargeToPeakVolts = 20;
1186 double ADCSignalConversionGain = ChargeToPeakVolts * 1.60e-04 * 2.4;
1187 double adc_input_voltage = input_electrons * ADCSignalConversionGain;
1188 unsigned int adc_output = (
unsigned int) (adc_input_voltage * 1024.0 / 2200.0);
1189 if (adc_output > 1023) adc_output = 1023;