7 #include <calobase/RawCluster.h>
8 #include <calobase/RawClusterContainer.h>
9 #include <calobase/RawTower.h>
10 #include <calobase/RawTowerContainer.h>
11 #include <calobase/RawTowerGeom.h>
12 #include <calobase/RawTowerGeomContainer.h>
24 #include <TLorentzVector.h>
26 #include <gsl/gsl_randist.h>
27 #include <gsl/gsl_rng.h>
34 return (a.second < b.second);
39 float deta = eta1 - eta2;
40 float dphi = phi1 - phi2;
41 while ( dphi > 3.14159 ) dphi -= 2 * 3.14159;
42 while ( dphi < -3.14159 ) dphi += 2 * 3.14159;
43 return sqrt( pow( deta, 2 ) + pow( dphi ,2 ) );
52 std::pair<float, float> expected_signature( response , resolution );
54 return expected_signature;
61 _energy_match_Nsigma( 1.5 ) ,
62 _emulate_efficiency( 1.1 )
65 _tr_eff = gsl_rng_alloc(gsl_rng_mt19937);
97 if (!pflowContainer) {
98 std::cout <<
" ERROR -- can't find ParticleFlowElements node after it should have been created" << std::endl;
103 int global_pflow_index = 0;
106 RawTowerContainer *towersEM = findNode::getClass<RawTowerContainer>(topNode,
"TOWER_CALIB_CEMC");
107 RawTowerContainer *towersIH = findNode::getClass<RawTowerContainer>(topNode,
"TOWER_CALIB_HCALIN");
108 RawTowerContainer *towersOH = findNode::getClass<RawTowerContainer>(topNode,
"TOWER_CALIB_HCALOUT");
110 if ( !towersEM || !towersIH || !towersOH ) {
111 std::cout <<
"ParticleFlowReco::process_event : FATAL ERROR, cannot find tower containers" << std::endl;
116 RawTowerGeomContainer *geomEM = findNode::getClass<RawTowerGeomContainer>(topNode,
"TOWERGEOM_CEMC");
117 RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode,
"TOWERGEOM_HCALIN");
118 RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode,
"TOWERGEOM_HCALOUT");
120 if ( !geomEM || !geomIH || !geomOH ) {
121 std::cout <<
"ParticleFlowReco::process_event : FATAL ERROR, cannot find tower geometry containers" << std::endl;
126 RawClusterContainer *clustersEM = findNode::getClass<RawClusterContainer>(topNode,
"TOPOCLUSTER_EMCAL");
127 RawClusterContainer *clustersHAD = findNode::getClass<RawClusterContainer>(topNode,
"TOPOCLUSTER_HCAL");
130 std::cout <<
"ParticleFlowReco::process_event : FATAL ERROR, cannot find cluster container TOPOCLUSTER_EMCAL" << std::endl;
133 if ( !clustersHAD ) {
134 std::cout <<
"ParticleFlowReco::process_event : FATAL ERROR, cannot find cluster container TOPOCLUSTER_HCAL" << std::endl;
164 std::cout <<
"ParticleFlowReco::process_event : initial population of TRK, EM, HAD objects " << std::endl;
179 float truth_pt = t.Pt();
180 float truth_p = t.P();
184 float truth_eta = t.Eta();
185 if (fabs (truth_eta) > 1.1)
188 float truth_phi = t.Phi();
193 if ( truth_pid != 211 && truth_pid != 321 && truth_pid != 2212 )
continue;
198 float this_eff = gsl_rng_uniform_pos(
_tr_eff);
217 std::cout <<
" TRK with p / pT = " << truth_p <<
" / " << truth_pt <<
" , eta / phi = " << truth_eta <<
" / " << truth_phi << std::endl;
230 float cluster_E = hiter->second->get_energy();
231 if ( cluster_E < 0.2 )
continue;
233 float cluster_phi = hiter->second->get_phi();
235 float cluster_theta = 3.14159 / 2.0 - atan2( hiter->second->get_z() , hiter->second->get_r() );
236 float cluster_eta = -1 * log( tan( cluster_theta / 2.0 ) );
245 if (
Verbosity() > 5 && cluster_E > 0.2 )
246 std::cout <<
" EM topoCluster with E = " << cluster_E <<
", eta / phi = " << cluster_eta <<
" / " << cluster_phi <<
" , nTow = " << hiter->second->getNTowers() << std::endl;
248 std::vector<float> this_cluster_tower_eta;
249 std::vector<float> this_cluster_tower_phi;
259 this_cluster_tower_phi.push_back( tower_geom->
get_phi() );
260 this_cluster_tower_eta.push_back( tower_geom->
get_eta() );
263 std::cout <<
"ParticleFlowReco::process_event : FATAL ERROR , EM topoClusters seem to contain HCal towers" << std::endl;
280 float cluster_E = hiter->second->get_energy();
281 if ( cluster_E < 0.2 )
continue;
283 float cluster_phi = hiter->second->get_phi();
285 float cluster_theta = 3.14159 / 2.0 - atan2( hiter->second->get_z() , hiter->second->get_r() );
286 float cluster_eta = -1 * log( tan( cluster_theta / 2.0 ) );
295 if (
Verbosity() > 5 && cluster_E > 0.2 )
296 std::cout <<
" HAD topoCluster with E = " << cluster_E <<
", eta / phi = " << cluster_eta <<
" / " << cluster_phi <<
" , nTow = " << hiter->second->getNTowers() << std::endl;
298 std::vector<float> this_cluster_tower_eta;
299 std::vector<float> this_cluster_tower_phi;
307 RawTower* tower = towersIH->getTower(iter->first);
310 this_cluster_tower_phi.push_back( tower_geom->
get_phi() );
311 this_cluster_tower_eta.push_back( tower_geom->
get_eta() );
316 RawTower* tower = towersOH->getTower(iter->first);
319 this_cluster_tower_phi.push_back( tower_geom->
get_phi() );
320 this_cluster_tower_eta.push_back( tower_geom->
get_eta() );
322 std::cout <<
"ParticleFlowReco::process_event : FATAL ERROR , HCal topoClusters seem to contain EM towers" << std::endl;
340 std::cout <<
"ParticleFlowReco::process_event : TRK -> EM and TRK -> HAD linking " << std::endl;
342 for (
unsigned int trk = 0; trk <
_pflow_TRK_p.size() ; trk++ ) {
348 float min_em_dR = 0.2;
349 int min_em_index = -1;
351 for (
unsigned int em = 0 ; em <
_pflow_EM_E.size() ; em++) {
354 if ( dR > 0.2 )
continue;
356 bool has_overlap =
false;
365 if ( dphi > 3.14159 ) dphi -= 2 * 3.14159;
366 if ( dphi < -3.14159 ) dphi += 2 * 3.14159;
368 if ( fabs( deta ) < 0.025 * 2.5 && fabs( dphi ) < 0.025 * 2.5 ) {
378 std::cout <<
" -> possible match to EM " << em <<
" with dR = " << dR << std::endl;
385 std::cout <<
" -> no match to EM " << em <<
" (even though dR = " << dR <<
" )" << std::endl;
408 if ( min_em_index > -1 ) {
413 std::cout <<
" -> matched EM " << min_em_index <<
" with pt / eta / phi = " <<
_pflow_EM_E.at( min_em_index ) <<
" / " <<
_pflow_EM_eta.at( min_em_index ) <<
" / " <<
_pflow_EM_phi.at( min_em_index ) <<
", dR = " << min_em_dR;
420 std::cout <<
" -> no EM match! ( best dR = " << min_em_dR <<
" ) " << std::endl;
424 float min_had_dR = 0.2;
425 int min_had_index = -1;
426 float max_had_pt = 0;
429 for (
unsigned int had = 0 ; had <
_pflow_HAD_E.size() ; had++) {
432 if ( dR > 0.5 )
continue;
434 bool has_overlap =
false;
443 if ( dphi > 3.14159 ) dphi -= 2 * 3.14159;
444 if ( dphi < -3.14159 ) dphi += 2 * 3.14159;
446 if ( fabs( deta ) < 0.1 * 1.5 && fabs( dphi ) < 0.1 * 1.5 ) {
456 std::cout <<
" -> possible match to HAD " << had <<
" with dR = " << dR << std::endl;
467 std::cout <<
" -> no match to HAD " << had <<
" (even though dR = " << dR <<
" )" << std::endl;
473 if ( min_had_index > -1 ) {
478 std::cout <<
" -> matched HAD " << min_had_index <<
" with pt / eta / phi = " <<
_pflow_HAD_E.at( min_had_index ) <<
" / " <<
_pflow_HAD_eta.at( min_had_index ) <<
" / " <<
_pflow_HAD_phi.at( min_had_index ) <<
", dR = " << min_had_dR << std::endl;
482 std::cout <<
" -> no HAD match! ( best dR = " << min_had_dR <<
" ) " << std::endl;
491 std::cout <<
"ParticleFlowReco::process_event : EM -> HAD linking " << std::endl;
493 for (
unsigned int em = 0; em <
_pflow_EM_E.size() ; em++ ) {
499 float min_had_dR = 0.2;
500 int min_had_index = -1;
501 float max_had_pt = 0;
503 for (
unsigned int had = 0 ; had <
_pflow_HAD_E.size() ; had++) {
506 if ( dR > 0.5 )
continue;
508 bool has_overlap =
false;
517 if ( dphi > 3.14159 ) dphi -= 2 * 3.14159;
518 if ( dphi < -3.14159 ) dphi += 2 * 3.14159;
520 if ( fabs( deta ) < 0.1 * 1.5 && fabs( dphi ) < 0.1 * 1.5 ) {
530 std::cout <<
" -> possible match to HAD " << had <<
" with dR = " << dR << std::endl;
541 std::cout <<
" -> no match to HAD " << had <<
" (even though dR = " << dR <<
" )" << std::endl;
547 if ( min_had_index > -1 ) {
552 std::cout <<
" -> matched HAD with E / eta / phi = " <<
_pflow_HAD_E.at( min_had_index ) <<
" / " <<
_pflow_HAD_eta.at( min_had_index ) <<
" / " <<
_pflow_HAD_phi.at( min_had_index ) <<
", dR = " << min_had_dR << std::endl;
556 std::cout <<
" -> no HAD match! ( best dR = " << min_had_dR <<
" ) " << std::endl;
564 std::cout <<
"ParticleFlowReco::process_event : sequential TRK -> EM && EM -> HAD ==> TRK -> HAD matching " << std::endl;
566 for (
unsigned int trk = 0; trk <
_pflow_TRK_p.size() ; trk++ ) {
579 bool is_trk_matched_to_HAD =
false;
582 if ( had == existing_had ) is_trk_matched_to_HAD =
true;
586 if ( ! is_trk_matched_to_HAD ) {
592 std::cout <<
" -> sequential match to HAD " << had <<
" through EM " << j << std::endl;
605 std::cout <<
"ParticleFlowReco::process_event : resolve TRK(s) + EM(s) -> HAD systems " << std::endl;
607 for (
unsigned int had = 0; had <
_pflow_HAD_E.size() ; had++ ) {
617 float total_TRK_p = 0;
618 float total_expected_E = 0;
619 float total_expected_E_var = 0;
654 float expected_E_mean = expected_signature.first;
655 float expected_E_sigma = expected_signature.second;
658 std::cout <<
" -> -> -> expected calo signature is " << expected_E_mean <<
" +/- " << expected_E_sigma << std::endl;
661 total_expected_E += expected_E_mean;
662 total_expected_E_var += pow( expected_E_sigma , 2 );
670 pflow->
set_px( tlv.Px() );
671 pflow->
set_py( tlv.Py() );
672 pflow->
set_pz( tlv.Pz() );
673 pflow->
set_e( tlv.E() );
674 pflow->
set_id( global_pflow_index );
675 pflow->
set_type( ParticleFlowElement::PFLOWTYPE::MATCHED_CHARGED_HADRON );
678 global_pflow_index++;
684 float total_expected_E_err = sqrt( total_expected_E_var );
687 std::cout <<
" -> Total track Sum p = " << total_TRK_p <<
" , expected calo Sum E = " << total_expected_E <<
" +/- " << total_expected_E_err <<
" , observed EM+HAD Sum E = " << total_EMHAD_E << std::endl;
692 if ( total_expected_E > total_EMHAD_E ) {
695 std::cout <<
" -> Expected E > Observed E, looking for additional potential TRK->EM matches" << std::endl;
697 std::map<int, float> additional_EMs;
706 std::cout <<
" -> -> TRK " << trk <<
" has " << addtl_matches <<
" additional matches! " << std::endl;
712 float existing_dR = 0.21;
726 std::vector< std::pair<int,float> > additional_EMs_vec;
728 for (
auto&
x : additional_EMs) {
729 additional_EMs_vec.push_back( std::pair<int,float>(
x.first ,
x.second ) );
735 std::cout <<
" -> Sorting the set of potential additional EMs " << std::endl;
740 while ( additional_EMs_vec.size() != 0 && total_expected_E > total_EMHAD_E ) {
742 int new_EM = additional_EMs_vec.at( 0 ).first;
745 std::cout <<
" -> adding EM " << new_EM <<
" ( dR = " << additional_EMs_vec.at( 0 ).second <<
" to the system (should not see it as orphan below)" << std::endl;
755 additional_EMs_vec.erase( additional_EMs_vec.begin() );
761 if ( n_EM_added > 0 ) {
762 std::cout <<
"After adding N = " << n_EM_added <<
" any additional EMs : " << std::endl;
763 std::cout <<
"-> Total track Sum p = " << total_TRK_p <<
" , expected calo Sum E = " << total_expected_E <<
" +/- " << total_expected_E_err <<
" , observed EM+HAD Sum E = " << total_EMHAD_E << std::endl;
766 std::cout <<
"No additional EMs found, continuing hypothesis check" << std::endl;
776 std::cout <<
" -> -> calo compatible within Nsigma = " <<
_energy_match_Nsigma <<
" , remove and keep tracks " << std::endl;
783 float residual_energy = total_EMHAD_E - total_expected_E;
786 std::cout <<
" -> -> calo not compatible, create leftover cluster with " << residual_energy << std::endl;
795 pflow->
set_px( tlv.Px() );
796 pflow->
set_py( tlv.Py() );
797 pflow->
set_pz( tlv.Pz() );
798 pflow->
set_e( tlv.E() );
799 pflow->
set_id( global_pflow_index );
800 pflow->
set_type( ParticleFlowElement::PFLOWTYPE::LEFTOVER_EM_PARTICLE );
803 global_pflow_index++;
812 std::cout <<
"ParticleFlowReco::process_event : resolve TRK(s) -> EM(s) ( + no HAD) systems " << std::endl;
814 for (
unsigned int em = 0; em <
_pflow_EM_E.size() ; em++ ) {
825 float total_TRK_p = 0;
826 float total_expected_E = 0;
827 float total_expected_E_var = 0;
845 float expected_E_mean = expected_signature.first;
846 float expected_E_sigma = expected_signature.second;
849 std::cout <<
" -> -> -> expected calo signature is " << expected_E_mean <<
" +/- " << expected_E_sigma << std::endl;
852 total_expected_E += expected_E_mean;
853 total_expected_E_var += pow( expected_E_sigma , 2 );
861 pflow->
set_px( tlv.Px() );
862 pflow->
set_py( tlv.Py() );
863 pflow->
set_pz( tlv.Pz() );
864 pflow->
set_e( tlv.E() );
865 pflow->
set_id( global_pflow_index );
866 pflow->
set_type( ParticleFlowElement::PFLOWTYPE::MATCHED_CHARGED_HADRON );
869 global_pflow_index++;
874 float total_expected_E_err = sqrt( total_expected_E_var );
877 std::cout <<
" -> Total track Sum p = " << total_TRK_p <<
" , expected calo Sum E = " << total_expected_E <<
" +/- " << total_expected_E_err <<
" , observed EM Sum E = " << total_EM_E << std::endl;
883 std::cout <<
" -> -> calo compatible within Nsigma = " <<
_energy_match_Nsigma <<
" , remove and keep tracks " << std::endl;
890 float residual_energy = total_EM_E - total_expected_E;
893 std::cout <<
" -> -> calo not compatible, create leftover cluster with " << residual_energy << std::endl;
902 pflow->
set_px( tlv.Px() );
903 pflow->
set_py( tlv.Py() );
904 pflow->
set_pz( tlv.Pz() );
905 pflow->
set_e( tlv.E() );
906 pflow->
set_id( global_pflow_index );
907 pflow->
set_type( ParticleFlowElement::PFLOWTYPE::LEFTOVER_EM_PARTICLE );
910 global_pflow_index++;
919 std::cout <<
"ParticleFlowReco::process_event : remove TRK-unlinked EMs and HADs " << std::endl;
921 for (
unsigned int em = 0; em <
_pflow_EM_E.size() ; em++ ) {
936 pflow->
set_px( tlv.Px() );
937 pflow->
set_py( tlv.Py() );
938 pflow->
set_pz( tlv.Pz() );
939 pflow->
set_e( tlv.E() );
940 pflow->
set_id( global_pflow_index );
941 pflow->
set_type( ParticleFlowElement::PFLOWTYPE::UNMATCHED_EM_PARTICLE );
944 global_pflow_index++;
949 for (
unsigned int had = 0; had <
_pflow_HAD_E.size() ; had++ ) {
964 pflow->
set_px( tlv.Px() );
965 pflow->
set_py( tlv.Py() );
966 pflow->
set_pz( tlv.Pz() );
967 pflow->
set_e( tlv.E() );
968 pflow->
set_id( global_pflow_index );
969 pflow->
set_type( ParticleFlowElement::PFLOWTYPE::UNMATCHED_NEUTRAL_HADRON );
972 global_pflow_index++;
977 for (
unsigned int trk = 0; trk <
_pflow_TRK_p.size() ; trk++ ) {
992 pflow->
set_px( tlv.Px() );
993 pflow->
set_py( tlv.Py() );
994 pflow->
set_pz( tlv.Pz() );
995 pflow->
set_e( tlv.E() );
996 pflow->
set_id( global_pflow_index );
997 pflow->
set_type( ParticleFlowElement::PFLOWTYPE::UNMATCHED_CHARGED_HADRON );
1000 global_pflow_index++;
1006 std::cout <<
"ParticleFlowReco::process_event: summary of PFlow objects " << std::endl;
1010 hiter->second->identify();
1026 std::cout <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
1035 dstNode->
addNode( pflowNode );
1039 ParticleFlowElementContainer *pflowElementContainer = findNode::getClass<ParticleFlowElementContainer>(topNode,
"ParticleFlowElements");
1040 if (!pflowElementContainer)
1044 pflowNode->
addNode( pflowElementNode );
1048 std::cout <<
PHWHERE <<
"::ERROR - ParticleFlowElements node alerady exists, but should not" << std::endl;
1060 std::cout <<
"ParticleFlowReco::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
1070 std::cout <<
"ParticleFlowReco::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
1080 std::cout <<
"ParticleFlowReco::End(PHCompositeNode *topNode) This is the End..." << std::endl;
1090 std::cout <<
"ParticleFlowReco::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
1098 std::cout <<
"ParticleFlowReco::Print(const std::string &what) const Printing info for " << what << std::endl;