16 #include <phparameter/PHParameterInterface.h>
40 #include <boost/foreach.hpp>
57 chkenergyconservation(0),
60 light_collection_model()
82 std::cout <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
89 cout <<
"PHG4FullProjSpacalCellReco::InitRun - Could not locate g4 hit node "
112 findNode::getClass<PHG4CylinderGeomContainer>(topNode,
geonodename.c_str());
115 cout <<
"PHG4FullProjSpacalCellReco::InitRun - Could not locate geometry node "
122 cout <<
"PHG4FullProjSpacalCellReco::InitRun - incoming geometry:"
137 assert(layergeom_raw);
145 cout <<
"PHG4FullProjSpacalCellReco::InitRun - Fatal Error -"
146 <<
" require to work with a version of SPACAL geometry of PHG4CylinderGeom_Spacalv3 or higher. "
147 <<
"However the incoming geometry has version "
148 << layergeom_raw->ClassName() << endl;
177 const double deltaphi = 2. *
M_PI / nphibin;
179 typedef map<double, int> map_z_tower_z_ID_t;
180 map_z_tower_z_ID_t map_z_tower_z_ID;
181 double phi_min = NAN;
183 BOOST_FOREACH(
const PHG4CylinderGeom_Spacalv3::tower_map_t::value_type& tower_pair, tower_map)
186 const int & tower_ID = tower_pair.first;
193 const int & tower_ID_z = tower_z_phi_ID.first;
194 const int & tower_ID_phi = tower_z_phi_ID.second;
196 if (tower_ID_phi == 0)
200 + sector_map.begin()->second;
206 map_z_tower_z_ID[tower.
centralZ] = tower_ID_z;
211 assert(! std::isnan(phi_min));
218 BOOST_FOREACH( map_z_tower_z_ID_t::value_type& z_tower_z_ID, map_z_tower_z_ID)
220 tower_z_ID_eta_bin_map[z_tower_z_ID.second] = eta_bin;
229 BOOST_FOREACH(
const PHG4CylinderGeom_Spacalv3::tower_map_t::value_type& tower_pair, tower_map)
232 const int & tower_ID = tower_pair.first;
238 const int & tower_ID_z = tower_z_phi_ID.first;
239 const int & tower_ID_phi = tower_z_phi_ID.second;
240 const int &
etabin = tower_z_ID_eta_bin_map[tower_ID_z];
248 auto z_to_eta = [&tower_radial](
const double &
z){
return -log(tan(0.5 * atan2(tower_radial,
z)));};
250 const double eta_central = z_to_eta(tower.
centralZ);
252 const double deta = (z_to_eta( tower.
centralZ + dz) - z_to_eta( tower.
centralZ - dz))/2;
255 for (
int sub_tower_ID_y = 0; sub_tower_ID_y < tower.
NSubtowerY;
260 const int sub_tower_etabin = etabin * layergeom->
get_n_subtower_eta() + sub_tower_ID_y;
262 const pair<double, double>etabounds (eta_central - deta + sub_tower_ID_y * 2 * deta / tower.
NSubtowerY,
263 eta_central - deta + (sub_tower_ID_y + 1) * 2 * deta / tower.
NSubtowerY);
265 const pair<double, double>zbounds (tower.
centralZ - dz + sub_tower_ID_y * 2 * dz / tower.
NSubtowerY,
269 layerseggeo->
set_zbounds(sub_tower_etabin,zbounds);
273 cout <<
"PHG4FullProjSpacalCellReco::InitRun::" <<
Name()
274 <<
"\t tower_ID_z = " << tower_ID_z
275 <<
"\t tower_ID_phi = " << tower_ID_phi
276 <<
"\t sub_tower_ID_y = " << sub_tower_ID_y
277 <<
"\t sub_tower_etabin = " << sub_tower_etabin
279 <<
"\t tower_radial = " << tower_radial
280 <<
"\t eta_central = " << eta_central
281 <<
"\t deta = " << deta
282 <<
"\t etabounds = [" <<etabounds.first <<
", " << etabounds.second<<
"]"
283 <<
"\t zbounds = [" <<zbounds.first <<
", " << zbounds.second<<
"]"
297 cout <<
"PHG4FullProjSpacalCellReco::InitRun::" <<
Name()
298 <<
" - Done layer" << (layergeom->
get_layer())
299 <<
". Print out the cell geometry:" << endl;
313 string paramnodename =
"G4CELLPARAM_" +
detector;
337 cout <<
"PHG4FullProjSpacalCellReco::process_event - Fatal Error - Could not locate g4 hit node "
344 cout <<
"PHG4FullProjSpacalCellReco::process_event - Fatal Error - could not locate cell node "
352 cout <<
"PHG4FullProjSpacalCellReco::process_event - Fatal Error - could not locate cell geometry node "
360 cout <<
"PHG4FullProjSpacalCellReco::process_event - Fatal Error - Could not locate sim geometry node "
373 assert(layergeom_raw);
379 for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
382 if (hiter->second->get_t(0)>
tmax)
continue;
383 if (hiter->second->get_t(1)<
tmin)
continue;
388 int scint_id = hiter->second->get_scint_id();
395 const int & tower_ID_z = tower_z_phi_ID.first;
396 const int & tower_ID_phi = tower_z_phi_ID.second;
398 PHG4CylinderGeom_Spacalv3::tower_map_t::const_iterator it_tower =
402 unsigned int key =
static_cast<unsigned int>(scint_id);
404 map<unsigned int, PHG4Cell *>::iterator
it =
celllist.find(key);
418 catch (exception &
e)
420 cout <<
"Print cell geometry:" << endl;
422 cout <<
"Print scint_id_coder:" << endl;
424 cout <<
"Print the hit:" << endl;
425 hiter->second->print();
426 cout <<
"PHG4FullProjSpacalCellReco::process_event::"
427 <<
Name() <<
" - Fatal Error - " << e.what() << endl;
431 const int sub_tower_ID_x = it_tower->second.get_sub_tower_ID_x(decoder.
fiber_ID);
432 const int sub_tower_ID_y = it_tower->second.get_sub_tower_ID_y(decoder.
fiber_ID);
433 unsigned short fiber_ID = decoder.
fiber_ID;
434 unsigned short etabinshort = etabin * layergeom->
get_n_subtower_eta() + sub_tower_ID_y;
435 unsigned short phibin = tower_ID_phi * layergeom->
get_n_subtower_phi() + sub_tower_ID_x;
447 * (hiter->second->get_local_z(0)
448 + hiter->second->get_local_z(1));
449 assert(not std::isnan(z));
457 const double x = it_tower->second.get_position_fraction_x_in_sub_tower(decoder.
fiber_ID);
458 const double y = it_tower->second.get_position_fraction_y_in_sub_tower(decoder.
fiber_ID);
463 cell->
add_edep(hiter->first, hiter->second->get_edep());
464 cell->
add_edep(hiter->second->get_edep());
466 cell->
add_shower_edep(hiter->second->get_shower_id(), hiter->second->get_edep());
470 for (map<unsigned int, PHG4Cell *>::const_iterator mapiter =
473 cells->
AddCell(mapiter->second);
477 cout <<
"PHG4FullProjSpacalCellReco::process_event::" <<
Name()
478 <<
" - " <<
"Adding cell in bin eta "
485 << mapiter->second->get_edep() <<
", light yield: "
486 << mapiter->second->get_light_yield() << endl;
492 cout <<
"PHG4FullProjSpacalCellReco::process_event::" <<
Name()
493 <<
" - " <<
" found " << numcells
494 <<
" fibers with energy deposition" << endl;
509 double sum_energy_cells = 0.;
512 for (citer = cell_begin_end.first; citer != cell_begin_end.second; ++citer)
514 sum_energy_cells += citer->second->get_edep();
521 <<
"PHG4FullProjSpacalCellReco::CheckEnergy - energy mismatch between cells: "
523 <<
" diff sum(cells) - sum(hits): "
531 cout <<
"PHG4FullProjSpacalCellReco::CheckEnergy::" <<
Name()
533 <<
" GeV. Passed CheckEnergy" << endl;
540 data_grid_light_guide_efficiency(nullptr),
541 data_grid_fiber_trans(nullptr)
545 "light collection efficiency as used in PHG4FullProjSpacalCellReco::LightCollectionModel;x positio fraction;y position fraction",
546 100, 0., 1., 100, 0., 1.);
549 "SCSF-78 Fiber Transmission as used in PHG4FullProjSpacalCellReco::LightCollectionModel;position in fiber (cm);Effective transmission",
562 delete data_grid_light_guide_efficiency;
563 delete data_grid_fiber_trans;
568 const std::string & input_file,
569 const std::string & histogram_light_guide_model,
570 const std::string & histogram_fiber_model)
572 TFile *
fin = TFile::Open(input_file.c_str());
575 assert(fin->IsOpen());
577 if (data_grid_light_guide_efficiency)
delete data_grid_light_guide_efficiency;
578 data_grid_light_guide_efficiency =
dynamic_cast<TH2 *
>(fin->Get(histogram_light_guide_model.c_str()));
579 assert(data_grid_light_guide_efficiency);
580 data_grid_light_guide_efficiency->SetDirectory(
nullptr);
582 if (data_grid_fiber_trans)
delete data_grid_fiber_trans;
583 data_grid_fiber_trans =
dynamic_cast<TH1 *
>(fin->Get(histogram_fiber_model.c_str()));
584 assert(data_grid_fiber_trans);
585 data_grid_fiber_trans->SetDirectory(
nullptr);
593 assert(data_grid_light_guide_efficiency);
594 assert(x_fraction >= 0);
595 assert(x_fraction <= 1);
596 assert(y_fraction >= 0);
597 assert(y_fraction <= 1);
599 const double eff = data_grid_light_guide_efficiency->Interpolate(x_fraction,
602 data_grid_light_guide_efficiency_verify->SetBinContent(
603 data_grid_light_guide_efficiency_verify->GetXaxis()->FindBin(x_fraction),
604 data_grid_light_guide_efficiency_verify->GetYaxis()->FindBin(y_fraction),
614 assert(data_grid_fiber_trans);
616 const double eff = data_grid_fiber_trans->Interpolate(z_distance);
618 data_grid_fiber_trans_verify->SetBinContent(
619 data_grid_fiber_trans_verify->GetXaxis()->FindBin(z_distance),