6 #include <calobase/RawTower.h>
7 #include <calobase/RawTowerContainer.h>
8 #include <calobase/RawTowerGeom.h>
9 #include <calobase/RawTowerGeomContainer.h>
28 #include <TLorentzVector.h>
49 , _backgroundName(
"TestTowerBackground")
54 _UE.resize(3, std::vector<float>(1, 0));
66 std::cout <<
"DetermineTowerBackground::process_event: entering with do_flow = " <<
_do_flow <<
", seed type = " <<
_seed_type <<
", ";
69 else if (_seed_type == 1)
72 std::cout <<
" UNDEFINED seed behavior! " << std::endl;
80 RawTowerContainer *towersEM3 = findNode::getClass<RawTowerContainer>(topNode,
"TOWER_CALIB_CEMC_RETOWER");
81 RawTowerContainer *towersIH3 = findNode::getClass<RawTowerContainer>(topNode,
"TOWER_CALIB_HCALIN");
82 RawTowerContainer *towersOH3 = findNode::getClass<RawTowerContainer>(topNode,
"TOWER_CALIB_HCALOUT");
85 std::cout <<
"DetermineTowerBackground::process_event: " << towersEM3->
size() <<
" TOWER_CALIB_CEMC_RETOWER towers" << std::endl;
86 std::cout <<
"DetermineTowerBackground::process_event: " << towersIH3->size() <<
" TOWER_CALIB_HCALIN towers" << std::endl;
87 std::cout <<
"DetermineTowerBackground::process_event: " << towersOH3->size() <<
" TOWER_CALIB_HCALOUT towers" << std::endl;
90 RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode,
"TOWERGEOM_HCALIN");
91 RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode,
"TOWERGEOM_HCALOUT");
96 JetMap *reco2_jets = findNode::getClass<JetMap>(topNode,
"AntiKt_Tower_HIRecoSeedsRaw_r02");
99 std::cout <<
"DetermineTowerBackground::process_event: examining possible seeds (1st iteration) ... " << std::endl;
103 Jet *this_jet = iter->second;
105 float this_pt = this_jet->
get_pt();
106 float this_phi = this_jet->
get_phi();
107 float this_eta = this_jet->
get_eta();
109 if (this_jet->
get_pt() < 5)
119 std::cout <<
"DetermineTowerBackground::process_event: possible seed jet with pt / eta / phi = " << this_pt <<
" / " << this_eta <<
" / " << this_phi <<
", examining constituents..." << std::endl;
121 std::map<int, double> constituent_ETsum;
132 if ((*comp).first == 5)
134 tower = towersIH3->getTower((*comp).second);
141 else if ((*comp).first == 7)
143 tower = towersOH3->getTower((*comp).second);
144 tower_geom = geomOH->get_tower_geometry(tower->
get_key());
150 else if ((*comp).first == 13)
152 tower = towersEM3->
getTower((*comp).second);
160 int comp_ikey = 1000 * comp_ieta + comp_iphi;
163 std::cout <<
"DetermineTowerBackground::process_event: --> --> constituent in layer " << (*comp).first <<
" at ieta / iphi = " << comp_ieta <<
" / " << comp_iphi <<
", filling map with key = " << comp_ikey <<
" and ET = " << comp_ET << std::endl;
165 constituent_ETsum[comp_ikey] += comp_ET;
168 std::cout <<
"DetermineTowerBackground::process_event: --> --> ET sum map at key = " << comp_ikey <<
" now has ET = " << constituent_ETsum[comp_ikey] << std::endl;
172 float constituent_max_ET = 0;
173 float constituent_sum_ET = 0;
174 int nconstituents = 0;
177 std::cout <<
"DetermineTowerBackground::process_event: --> now iterating over map..." << std::endl;
178 for (std::map<int, double>::iterator map_iter = constituent_ETsum.begin(); map_iter != constituent_ETsum.end(); ++map_iter)
181 std::cout <<
"DetermineTowerBackground::process_event: --> --> map has key # " << map_iter->first <<
" and ET = " << map_iter->second << std::endl;
183 constituent_sum_ET += map_iter->second;
184 if (map_iter->second > constituent_max_ET) constituent_max_ET = map_iter->second;
187 float mean_constituent_ET = constituent_sum_ET / nconstituents;
188 float seed_D = constituent_max_ET / mean_constituent_ET;
191 this_jet->
set_property(Jet::PROPERTY::prop_SeedD, seed_D);
194 std::cout <<
"DetermineTowerBackground::process_event: --> jet has < ET > = " << constituent_sum_ET <<
" / " << nconstituents <<
" = " << mean_constituent_ET <<
", max-ET = " << constituent_max_ET <<
", and D = " << seed_D << std::endl;
202 this_jet->
set_property(Jet::PROPERTY::prop_SeedItr, 1.0);
205 std::cout <<
"DetermineTowerBackground::process_event: --> adding seed at eta / phi = " << this_eta <<
" / " << this_phi <<
" ( R=0.2 jet with pt = " << this_pt <<
", D = " << seed_D <<
" ) " << std::endl;
210 this_jet->
set_property(Jet::PROPERTY::prop_SeedItr, 0.0);
213 std::cout <<
"DetermineTowerBackground::process_event: --> discarding potential seed at eta / phi = " << this_eta <<
" / " << this_phi <<
" ( R=0.2 jet with pt = " << this_pt <<
", D = " << seed_D <<
" ) " << std::endl;
223 JetMap *reco2_jets = findNode::getClass<JetMap>(topNode,
"AntiKt_Tower_HIRecoSeedsSub_r02");
226 std::cout <<
"DetermineTowerBackground::process_event: examining possible seeds (2nd iteration) ... " << std::endl;
230 Jet *this_jet = iter->second;
232 float this_pt = this_jet->
get_pt();
233 float this_phi = this_jet->
get_phi();
234 float this_eta = this_jet->
get_eta();
239 this_jet->
set_property(Jet::PROPERTY::prop_SeedItr, 0.0);
248 this_jet->
set_property(Jet::PROPERTY::prop_SeedItr, 2.0);
251 std::cout <<
"DetermineTowerBackground::process_event: --> adding seed at eta / phi = " << this_eta <<
" / " << this_phi <<
" ( R=0.2 jet with pt = " << this_pt <<
" ) " << std::endl;
278 std::cout <<
"DetermineTowerBackground::process_event: setting number of towers in eta / phi: " <<
_HCAL_NETA <<
" / " <<
_HCAL_NPHI << std::endl;
307 float this_eta = tower_geom->
get_eta();
308 float this_phi = tower_geom->
get_phi();
309 int this_etabin = geomIH->
get_etabin(this_eta);
310 int this_phibin = geomIH->
get_phibin(this_phi);
313 _EMCAL_E[this_etabin][this_phibin] += this_E;
317 std::cout <<
"DetermineTowerBackground::process_event: EMCal tower eta ( bin ) / phi ( bin ) / E = " << std::setprecision(6) << this_eta <<
" ( " << this_etabin <<
" ) / " << this_phi <<
" ( " << this_phibin <<
" ) / " << this_E << std::endl;
328 float this_eta = tower_geom->
get_eta();
329 float this_phi = tower_geom->
get_phi();
330 int this_etabin = geomIH->
get_etabin(this_eta);
331 int this_phibin = geomIH->
get_phibin(this_phi);
334 _IHCAL_E[this_etabin][this_phibin] += this_E;
338 std::cout <<
"DetermineTowerBackground::process_event: IHCal tower at eta ( bin ) / phi ( bin ) / E = " << std::setprecision(6) << this_eta <<
" ( " << this_etabin <<
" ) / " << this_phi <<
" ( " << this_phibin <<
" ) / " << this_E << std::endl;
349 float this_eta = tower_geom->
get_eta();
350 float this_phi = tower_geom->
get_phi();
351 int this_etabin = geomOH->get_etabin(this_eta);
352 int this_phibin = geomOH->get_phibin(this_phi);
355 _OHCAL_E[this_etabin][this_phibin] += this_E;
359 std::cout <<
"DetermineTowerBackground::process_event: OHCal tower at eta ( bin ) / phi ( bin ) / E = " << std::setprecision(6) << this_eta <<
" ( " << this_etabin <<
" ) / " << this_phi <<
" ( " << this_phibin <<
" ) / " << this_E << std::endl;
373 std::cout <<
"DetermineTowerBackground::process_event: flow not enabled, setting Psi2 = " <<
_Psi2 <<
" ( " <<
_Psi2 / 3.14159 <<
" * pi ) , v2 = " <<
_v2 << std::endl;
380 int nStripsAvailableForFlow = 0;
381 int nStripsUnavailableForFlow = 0;
388 for (
int eta = 0;
eta < local_max_eta;
eta++)
390 bool isAnyTowerExcluded =
false;
392 for (
int phi = 0;
phi < local_max_phi;
phi++)
397 bool isExcluded =
false;
399 for (
unsigned int iseed = 0; iseed <
_seed_eta.size(); iseed++)
401 float deta = this_eta -
_seed_eta[iseed];
402 float dphi = this_phi -
_seed_phi[iseed];
403 if (dphi > 3.14159) dphi -= 2 * 3.14159;
404 if (dphi < -3.14159) dphi += 2 * 3.14159;
405 float dR = sqrt(pow(deta, 2) + pow(dphi, 2));
413 std::cout <<
" setting excluded mark for tower with E / eta / phi = " << my_E <<
" / " << this_eta <<
" / " << this_phi <<
" from seed at eta / phi = " << _seed_eta[iseed] <<
" / " << _seed_phi[iseed] << std::endl;
421 isAnyTowerExcluded =
true;
425 if (!isAnyTowerExcluded)
428 std::cout <<
" strip at layer " <<
layer <<
", eta " <<
eta <<
" has no excluded towers and can be used for flow determination " << std::endl;
429 nStripsAvailableForFlow++;
431 for (
int phi = 0;
phi < local_max_phi;
phi++)
445 std::cout <<
" strip at layer " <<
layer <<
", eta " <<
eta <<
" DOES have excluded towers and CANNOT be used for flow determination " << std::endl;
446 nStripsUnavailableForFlow++;
459 std::cout <<
"DetermineTowerBackground::process_event: # of strips (summed over layers) available / unavailable for flow determination: " << nStripsAvailableForFlow <<
" / " << nStripsUnavailableForFlow << std::endl;
461 if (nStripsAvailableForFlow > 0)
472 _Psi2 = atan2(Q_y, Q_x) / 2.0;
480 std::cout <<
"DetermineTowerBackground::process_event: FATAL , G4TruthInfo does not exist , cannot extract truth flow with do_flow = " <<
_do_flow << std::endl;
486 float Hijing_Qx = 0, Hijing_Qy = 0;
497 float truth_pt = t.Pt();
498 if (truth_pt < 0.4)
continue;
499 float truth_eta = t.Eta();
500 if (fabs(truth_eta) > 1.1)
continue;
501 float truth_phi = t.Phi();
502 int truth_pid = g4particle->
get_pid();
505 std::cout <<
"DetermineTowerBackground::process_event: determining truth flow, using particle w/ pt / eta / phi " << truth_pt <<
" / " << truth_eta <<
" / " << truth_phi <<
" , embed / PID = " << truthinfo->
isEmbeded(g4particle->
get_track_id()) <<
" / " << truth_pid << std::endl;
507 Hijing_Qx += truth_pt * cos(2 * truth_phi);
508 Hijing_Qy += truth_pt * sin(2 * truth_phi);
511 _Psi2 = atan2(Hijing_Qy, Hijing_Qx) / 2.0;
514 std::cout <<
"DetermineTowerBackground::process_event: flow extracted from Hijing truth particles, setting Psi2 = " <<
_Psi2 <<
" ( " <<
_Psi2 / 3.14159 <<
" * pi ) " << std::endl;
518 double sum_cos2dphi = 0;
524 _v2 = sum_cos2dphi /
E;
534 std::cout <<
"DetermineTowerBackground::process_event: no full strips available for flow modulation, setting v2 and Psi = 0" << std::endl;
539 std::cout <<
"DetermineTowerBackground::process_event: unnormalized Q vector (Qx, Qy) = ( " << Q_x <<
", " << Q_y <<
" ) with Sum E_i = " << E << std::endl;
540 std::cout <<
"DetermineTowerBackground::process_event: Psi2 = " <<
_Psi2 <<
" ( " <<
_Psi2 / 3.14159 <<
" * pi " << (
_do_flow == 2 ?
"from Hijing " :
"") <<
") , v2 = " <<
_v2 <<
" ( using " <<
_nStrips <<
" ) " << std::endl;
553 for (
int eta = 0;
eta < local_max_eta;
eta++)
558 for (
int phi = 0;
phi < local_max_phi;
phi++)
563 bool isExcluded =
false;
565 for (
unsigned int iseed = 0; iseed <
_seed_eta.size(); iseed++)
567 float deta = this_eta -
_seed_eta[iseed];
568 float dphi = this_phi -
_seed_phi[iseed];
569 if (dphi > 3.14159) dphi -= 2 * 3.14159;
570 if (dphi < -3.14159) dphi += 2 * 3.14159;
571 float dR = sqrt(pow(deta, 2) + pow(dphi, 2));
575 if (
Verbosity() > 10) std::cout <<
" setting excluded mark from seed at eta / phi = " << _seed_eta[iseed] <<
" / " << _seed_phi[iseed] << std::endl;
589 if (
Verbosity() > 10) std::cout <<
" tower at eta / phi = " << this_eta <<
" / " << this_phi <<
" with E = " << total_E <<
" excluded due to seed " << std::endl;
594 std::pair<float, float> phibounds = geomIH->
get_phibounds(0);
596 float deta = etabounds.second - etabounds.first;
597 float dphi = phibounds.second - phibounds.first;
598 float total_area = total_tower * deta * dphi;
603 std::cout <<
"DetermineTowerBackground::process_event: at layer / eta index ( eta range ) = " <<
layer <<
" / " <<
eta <<
" ( " << etabounds.first <<
" - " << etabounds.second <<
" ) , total E / total Ntower / total area = " << total_E <<
" / " << total_tower <<
" / " << total_area <<
" , UE per tower = " << total_E / total_tower << std::endl;
612 std::cout <<
"DetermineTowerBackground::process_event: summary UE in layer " <<
layer <<
" : ";
614 std::cout << std::endl;
622 if (
Verbosity() > 0) std::cout <<
"DetermineTowerBackground::process_event: exiting" << std::endl;
635 std::cout <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
649 if (!towerbackground)
667 if (!towerbackground)
669 std::cout <<
" ERROR -- can't find TowerBackground node after it should have been created" << std::endl;