packages/jetbackground/DetermineTowerBackground.cc

0001 #include "DetermineTowerBackground.h"
0002
0003 #include "TowerBackground.h"
0004 #include "TowerBackgroundv1.h"
0005
0006 #include <calobase/RawTower.h>
0007 #include <calobase/RawTowerContainer.h>
0008 #include <calobase/RawTowerDefs.h>
0009 #include <calobase/RawTowerGeom.h>
0010 #include <calobase/RawTowerGeomContainer.h>
0011 #include <calobase/TowerInfo.h>
0012 #include <calobase/TowerInfoContainer.h>
0013
0014 #include <eventplaneinfo/Eventplaneinfo.h>
0015 #include <eventplaneinfo/EventplaneinfoMap.h>
0016
0017 #include <jetbase/Jet.h>
0018 #include <jetbase/JetContainer.h>
0019
0020 #include <g4main/PHG4Particle.h>
0021 #include <g4main/PHG4TruthInfoContainer.h>
0022
0023 #include <fun4all/Fun4AllReturnCodes.h>
0024 #include <fun4all/SubsysReco.h>
0025
0026 #include <phool/PHCompositeNode.h>
0027 #include <phool/PHIODataNode.h>
0028 #include <phool/PHNode.h>
0029 #include <phool/PHNodeIterator.h>
0030 #include <phool/PHObject.h>
0031 #include <phool/getClass.h>
0032 #include <phool/phool.h>
0033
0034 #include <TLorentzVector.h>
0035
0036 // standard includes
0037 #include <algorithm>
0038 #include <cmath>
0039 #include <cstdlib>
0040 #include <iomanip>
0041 #include <iostream>
0042 #include <map>
0043 #include <utility>
0044 #include <vector>
0045 #include <set>
0046
0047 DetermineTowerBackground::DetermineTowerBackground(const std::string &name)
0048   : SubsysReco(name)
0049 {
0050   _UE.resize(3, std::vector<float>(1, 0));
0051 }
0052
0053 int DetermineTowerBackground::InitRun(PHCompositeNode *topNode)
0054 {
0055   return CreateNode(topNode);
0056 }
0057
0058 int DetermineTowerBackground::process_event(PHCompositeNode *topNode)
0059 {
0060
0061
0062   if ( Verbosity() > 0 )
0063   {
0064     std::cout << "DetermineTowerBackground::process_event: entering with do_flow = " << _do_flow << ", seed type = " << _seed_type << ", ";
0065     if (_seed_type == 0)
0066     {
0067       std::cout << " D = " << _seed_jet_D << std::endl;
0068     }
0069     else if (_seed_type == 1)
0070     {
0071       std::cout << " pT = " << _seed_jet_pt << std::endl;
0072     }
0073     else
0074     {
0075       std::cout << " UNDEFINED seed behavior! " << std::endl;
0076     }
0077   }
0078
0079
0080   // pull out the tower containers and geometry objects at the start
0081   EMTowerName = m_towerNodePrefix + "_CEMC_RETOWER";
0082   IHTowerName = m_towerNodePrefix + "_HCALIN";
0083   OHTowerName = m_towerNodePrefix + "_HCALOUT";
0084   auto * towerinfosEM3 = findNode::getClass<TowerInfoContainer>(topNode, EMTowerName);
0085   auto * towerinfosIH3 = findNode::getClass<TowerInfoContainer>(topNode, IHTowerName);
0086   auto * towerinfosOH3 = findNode::getClass<TowerInfoContainer>(topNode, OHTowerName);
0087   if ( !towerinfosEM3 )
0088   {
0089     std::cout << "DetermineTowerBackground::process_event: Cannot find node " << EMTowerName << std::endl;
0090     exit(1);
0091   }
0092   if ( !towerinfosIH3 )
0093   {
0094     std::cout << "DetermineTowerBackground::process_event: Cannot find node " << IHTowerName << std::endl;
0095     exit(1);
0096   }
0097   if ( !towerinfosOH3 )
0098   {
0099     std::cout << "DetermineTowerBackground::process_event: Cannot find node " << OHTowerName << std::endl;
0100     exit(1);
0101   }
0102
0103   auto * geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0104   auto * geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0105   if ( !geomIH )
0106   {
0107     std::cout << "DetermineTowerBackground::process_event: Cannot find TOWERGEOM_HCALIN, exiting" << std::endl;
0108     exit(1);
0109   }
0110   if ( !geomOH )
0111   {
0112     std::cout << "DetermineTowerBackground::process_event: Cannot find TOWERGEOM_HCALOUT, exiting" << std::endl;
0113     exit(1);
0114   }
0115
0116   // clear seed eta/phi positions
0117   _seed_eta.resize(0);
0118   _seed_phi.resize(0);
0119
0120   // get the binning from the geometry (different for 1D vs 2D...)
0121   if ( _HCAL_NETA < 0 ) // fisrt event
0122   {
0123     _HCAL_NETA = geomIH->get_etabins();
0124     _HCAL_NPHI = geomIH->get_phibins();
0125
0126     // resize UE density and energy vectors
0127     _UE.resize(3 , std::vector<float>(_HCAL_NETA, 0));
0128
0129     _EMCAL_E.resize(_HCAL_NETA, std::vector<float>(_HCAL_NPHI, 0));
0130     _IHCAL_E.resize(_HCAL_NETA, std::vector<float>(_HCAL_NPHI, 0));
0131     _OHCAL_E.resize(_HCAL_NETA, std::vector<float>(_HCAL_NPHI, 0));
0132
0133     _EMCAL_ISBAD.resize(_HCAL_NETA, std::vector<int>(_HCAL_NPHI, 0));
0134     _IHCAL_ISBAD.resize(_HCAL_NETA, std::vector<int>(_HCAL_NPHI, 0));
0135     _OHCAL_ISBAD.resize(_HCAL_NETA, std::vector<int>(_HCAL_NPHI, 0));
0136
0137     // for flow determination, build up a 1-D phi distribution of
0138     // energies from all layers summed together, populated only from eta
0139     // strips which do not have any excluded phi towers
0140     _FULLCALOFLOW_PHI_E.resize(_HCAL_NPHI, 0);
0141     _FULLCALOFLOW_PHI_VAL.resize(_HCAL_NPHI, 0);
0142
0143     // defualt set weights to 1.0 for all phi bins
0144     _EMCAL_PHI_WEIGHTS.resize(_HCAL_NPHI, 1.0);
0145     _IHCAL_PHI_WEIGHTS.resize(_HCAL_NPHI, 1.0);
0146     _OHCAL_PHI_WEIGHTS.resize(_HCAL_NPHI, 1.0);
0147
0148     if (Verbosity() > 0)
0149     {
0150       std::cout << "DetermineTowerBackground::process_event: setting number of towers in eta / phi: " << _HCAL_NETA << " / " << _HCAL_NPHI << std::endl;
0151     }
0152   }
0153
0154   // reset all maps map
0155   _UE.assign(3, std::vector<float>(_HCAL_NETA, 0));
0156
0157   // reset all energy vectors
0158   _EMCAL_E.assign(_HCAL_NETA, std::vector<float>(_HCAL_NPHI, 0));
0159   _IHCAL_E.assign(_HCAL_NETA, std::vector<float>(_HCAL_NPHI, 0));
0160   _OHCAL_E.assign(_HCAL_NETA, std::vector<float>(_HCAL_NPHI, 0));
0161
0162   // reset bad tower masks
0163   _EMCAL_ISBAD.assign(_HCAL_NETA, std::vector<int>(_HCAL_NPHI, 0));
0164   _IHCAL_ISBAD.assign(_HCAL_NETA, std::vector<int>(_HCAL_NPHI, 0));
0165   _OHCAL_ISBAD.assign(_HCAL_NETA, std::vector<int>(_HCAL_NPHI, 0));
0166
0167   // create a set for all eta strips to be updated
0168   std::set<int> EtaStripsAvailbleForFlow = {};
0169   for ( int eta = 0; eta < _HCAL_NETA; eta++) { EtaStripsAvailbleForFlow.insert(eta); }
0170
0171   // seed type 0 is D > 3 R=0.2 jets run on retowerized CEMC
0172   if (_seed_type == 0)
0173   {
0174     auto * reco2_jets = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsRaw_r02");
0175     if (!reco2_jets)
0176     {
0177       std::cout << "DetermineTowerBackground::process_event: Cannot find AntiKt_TowerInfo_HIRecoSeedsRaw_r02, exiting" << std::endl;
0178       exit(1);
0179     }
0180     if (Verbosity() > 1)
0181     {
0182       std::cout << "DetermineTowerBackground::process_event: examining possible seeds (1st iteration) ... " << std::endl;
0183     }
0184
0185     _index_SeedD = reco2_jets->property_index(Jet::PROPERTY::prop_SeedD);
0186     _index_SeedItr = reco2_jets->property_index(Jet::PROPERTY::prop_SeedItr);
0187     for (auto *this_jet : *reco2_jets)
0188     {
0189       float this_pt = this_jet->get_pt();
0190       float this_phi = this_jet->get_phi();
0191       float this_eta = this_jet->get_eta();
0192
0193       if (this_jet->get_pt() < 5)
0194       {
0195         // mark that this jet was not selected as a seed (and did not have D determined)
0196         this_jet->set_property(_index_SeedD, 0);
0197         this_jet->set_property(_index_SeedItr, 0);
0198
0199         continue;
0200       }
0201
0202       if (Verbosity() > 2)
0203       {
0204         std::cout << "DetermineTowerBackground::process_event: possible seed jet with pt / eta / phi = " << this_pt << " / " << this_eta << " / " << this_phi << ", examining constituents..." << std::endl;
0205       }
0206
0207       std::map<int, double> constituent_ETsum;
0208
0209       for (const auto &comp : this_jet->get_comp_vec())
0210       {
0211         int comp_ieta = -1;
0212         int comp_iphi = -1;
0213         float comp_ET = 0;
0214         int comp_isBad = -99;
0215
0216         TowerInfo *towerinfo;
0217         RawTowerGeom *tower_geom;
0218
0219
0220         if (comp.first == 5 || comp.first == 26)
0221         {
0222           towerinfo = towerinfosIH3->get_tower_at_channel(comp.second);
0223           unsigned int towerkey = towerinfosIH3->encode_key(comp.second);
0224           comp_ieta = towerinfosIH3->getTowerEtaBin(towerkey);
0225           comp_iphi = towerinfosIH3->getTowerPhiBin(towerkey);
0226           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, comp_ieta, comp_iphi);
0227           tower_geom = geomIH->get_tower_geometry(key);
0228           comp_ET = towerinfo->get_energy() / cosh(tower_geom->get_eta());
0229           comp_isBad = towerinfo->get_isHot() || towerinfo->get_isNoCalib() || towerinfo->get_isNotInstr() || towerinfo->get_isBadChi2();
0230         }
0231         else if (comp.first == 7 || comp.first == 27)
0232         {
0233           towerinfo = towerinfosOH3->get_tower_at_channel(comp.second);
0234           unsigned int towerkey = towerinfosOH3->encode_key(comp.second);
0235           comp_ieta = towerinfosOH3->getTowerEtaBin(towerkey);
0236           comp_iphi = towerinfosOH3->getTowerPhiBin(towerkey);
0237           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, comp_ieta, comp_iphi);
0238           tower_geom = geomOH->get_tower_geometry(key);
0239           comp_ET = towerinfo->get_energy() / cosh(tower_geom->get_eta());
0240           comp_isBad = towerinfo->get_isHot() || towerinfo->get_isNoCalib() || towerinfo->get_isNotInstr() || towerinfo->get_isBadChi2();
0241         }
0242         else if (comp.first == 13 || comp.first == 28)
0243         {
0244           towerinfo = towerinfosEM3->get_tower_at_channel(comp.second);
0245           unsigned int towerkey = towerinfosEM3->encode_key(comp.second);
0246           comp_ieta = towerinfosEM3->getTowerEtaBin(towerkey);
0247           comp_iphi = towerinfosEM3->getTowerPhiBin(towerkey);
0248           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, comp_ieta, comp_iphi);
0249           tower_geom = geomIH->get_tower_geometry(key);
0250           comp_ET = towerinfo->get_energy() / cosh(tower_geom->get_eta());
0251           comp_isBad = towerinfo->get_isHot() || towerinfo->get_isNoCalib() || towerinfo->get_isNotInstr() || towerinfo->get_isBadChi2();
0252         }
0253
0254
0255         if (comp_isBad)
0256         {
0257           if (Verbosity() > 4)
0258           {
0259             std::cout << "DetermineTowerBackground::process_event: --> --> Skipping constituent in layer " << comp.first << " at ieta / iphi = " << comp_ieta << " / " << comp_iphi << "due to masking" << std::endl;
0260           }
0261           continue;
0262         }
0263         int comp_ikey = (1000 * comp_ieta) + comp_iphi;
0264
0265         if (Verbosity() > 4)
0266         {
0267           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;
0268         }
0269
0270         constituent_ETsum[comp_ikey] += comp_ET;
0271
0272         if (Verbosity() > 4)
0273         {
0274           std::cout << "DetermineTowerBackground::process_event: --> --> ET sum map at key = " << comp_ikey << " now has ET = " << constituent_ETsum[comp_ikey] << std::endl;
0275         }
0276       }
0277
0278       // now iterate over constituent_ET sums to find maximum and mean
0279       float constituent_max_ET = 0;
0280       float constituent_sum_ET = 0;
0281       int nconstituents = 0;
0282
0283       if (Verbosity() > 4)
0284       {
0285         std::cout << "DetermineTowerBackground::process_event: --> now iterating over map..." << std::endl;
0286       }
0287       for (auto &map_iter : constituent_ETsum)
0288       {
0289         if (Verbosity() > 4)
0290         {
0291           std::cout << "DetermineTowerBackground::process_event: --> --> map has key # " << map_iter.first << " and ET = " << map_iter.second << std::endl;
0292         }
0293         nconstituents++;
0294         constituent_sum_ET += map_iter.second;
0295         constituent_max_ET = std::max<double>(map_iter.second, constituent_max_ET);
0296       }
0297
0298       float mean_constituent_ET = constituent_sum_ET / nconstituents;
0299       float seed_D = constituent_max_ET / mean_constituent_ET;
0300
0301       // store D value as property for offline analysis / debugging
0302       this_jet->set_property(_index_SeedD, seed_D);
0303
0304       if (Verbosity() > 3)
0305       {
0306         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;
0307       }
0308
0309       if (seed_D > _seed_jet_D && constituent_max_ET > _seed_max_const)  // this will be constituent_max_ET > 0 if not set
0310       {
0311         _seed_eta.push_back(this_eta);
0312         _seed_phi.push_back(this_phi);
0313         int seed_ieta = geomIH->get_etabin(this_eta);
0314
0315         // remove eta-4 to eta+4 from the set of all eta strips
0316         // dont need to worry about bounds since itsa set
0317         for ( int ieta = -4; ieta <= 4; ieta++ ) { EtaStripsAvailbleForFlow.erase(seed_ieta + ieta); }
0318
0319         // set first iteration seed property
0320         this_jet->set_property(_index_SeedItr, 1.0);
0321
0322         if (Verbosity() > 1)
0323         {
0324           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;
0325         }
0326       }
0327       else
0328       {
0329         // mark that this jet was considered but not used as a seed
0330         this_jet->set_property(_index_SeedItr, 0.0);
0331
0332         if (Verbosity() > 3)
0333         {
0334           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;
0335         }
0336       }
0337     }
0338   }
0339
0340   // seed type 1 is the set of those jets above which, when their
0341   // kinematics are updated for the first background subtraction, have
0342   // pT > 20 GeV
0343   if (_seed_type == 1)
0344   {
0345     auto * reco2_jets = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsSub_r02");
0346     if (!reco2_jets)
0347     {
0348       std::cout << "DetermineTowerBackground::process_event: Cannot find AntiKt_TowerInfo_HIRecoSeedsSub_r02, exiting" << std::endl;
0349       exit(1);
0350     }
0351     if (Verbosity() > 1)
0352     {
0353       std::cout << "DetermineTowerBackground::process_event: examining possible seeds (2nd iteration) ... " << std::endl;
0354     }
0355
0356     _index_SeedD = reco2_jets->property_index(Jet::PROPERTY::prop_SeedD);
0357     _index_SeedItr = reco2_jets->property_index(Jet::PROPERTY::prop_SeedItr);
0358     for (auto *this_jet : *reco2_jets)
0359     {
0360       float this_pt = this_jet->get_pt();
0361       float this_phi = this_jet->get_phi();
0362       float this_eta = this_jet->get_eta();
0363
0364       if (this_jet->get_pt() < _seed_jet_pt)
0365       {
0366         // mark that this jet was considered but not used as a seed
0367         this_jet->set_property(_index_SeedItr, 0.0);
0368
0369         continue;
0370       }
0371
0372       _seed_eta.push_back(this_eta);
0373       _seed_phi.push_back(this_phi);
0374
0375       int seed_ieta = geomIH->get_etabin(this_eta);
0376
0377       // remove eta-4 to eta+4 from the set of all eta strips
0378       // dont need to worry about bounds since itsa set
0379       for ( int ieta = -4; ieta <= 4; ieta++ ) { EtaStripsAvailbleForFlow.erase(seed_ieta + ieta); }
0380
0381
0382       // set second iteration seed property
0383       this_jet->set_property(_index_SeedItr, 2.0);
0384
0385       if (Verbosity() > 1)
0386       {
0387         std::cout << "DetermineTowerBackground::process_event: --> adding seed at eta / phi = " << this_eta << " / " << this_phi << " ( R=0.2 jet with pt = " << this_pt << " ) " << std::endl;
0388       }
0389     }
0390   }
0391
0392
0393   int MaxEtaBinsWithoutSeeds = EtaStripsAvailbleForFlow.size();
0394   if (Verbosity() > 1)
0395   {
0396     for (const auto &eta : EtaStripsAvailbleForFlow)
0397     {
0398       std::cout << "DetermineTowerBackground::process_event: Remaining eta strip for background determination: " << eta << std::endl;
0399     }
0400     std::cout << "DetermineTowerBackground::process_event: Finished processing seeds. Remaining avilable eta strips for background determination: " << MaxEtaBinsWithoutSeeds << std::endl;
0401   }
0402
0403   // fill energy and status vectors
0404   // iterate over EMCal towerinfos
0405   if (!towerinfosEM3 || !towerinfosIH3 || !towerinfosOH3)
0406   {
0407     std::cout << PHWHERE << "missing tower info object, doing nothing" << std::endl;
0408     return Fun4AllReturnCodes::ABORTRUN;
0409   }
0410   unsigned int nchannels_em = towerinfosEM3->size();
0411   for (unsigned int channel = 0; channel < nchannels_em; channel++)
0412   {
0413     unsigned int key = towerinfosEM3->encode_key(channel);
0414     int this_etabin = towerinfosEM3->getTowerEtaBin(key);
0415     int this_phibin = towerinfosEM3->getTowerPhiBin(key);
0416     TowerInfo *tower = towerinfosEM3->get_tower_at_channel(channel);
0417     float this_E = tower->get_energy();
0418     int this_isBad = tower->get_isHot() || tower->get_isNoCalib() || tower->get_isNotInstr() || tower->get_isBadChi2();
0419     _EMCAL_ISBAD[this_etabin][this_phibin] = this_isBad;
0420     if (!this_isBad)
0421     { // just in case since all energy is summed
0422       _EMCAL_E[this_etabin][this_phibin] += this_E;
0423     }
0424
0425   }
0426
0427   // iterate over IHCal towerinfos
0428   unsigned int nchannels_ih = towerinfosIH3->size();
0429   for (unsigned int channel = 0; channel < nchannels_ih; channel++)
0430   {
0431     unsigned int key = towerinfosIH3->encode_key(channel);
0432     int this_etabin = towerinfosIH3->getTowerEtaBin(key);
0433     int this_phibin = towerinfosIH3->getTowerPhiBin(key);
0434     TowerInfo *tower = towerinfosIH3->get_tower_at_channel(channel);
0435     float this_E = tower->get_energy();
0436     int this_isBad = tower->get_isHot() || tower->get_isNoCalib() || tower->get_isNotInstr() || tower->get_isBadChi2();
0437     _IHCAL_ISBAD[this_etabin][this_phibin] = this_isBad;
0438     if (!this_isBad)
0439     { // just in case since all energy is summed
0440       _IHCAL_E[this_etabin][this_phibin] += this_E;
0441     }
0442
0443   }
0444
0445   // iterate over OHCal towerinfos
0446   unsigned int nchannels_oh = towerinfosOH3->size();
0447   for (unsigned int channel = 0; channel < nchannels_oh; channel++)
0448   {
0449     unsigned int key = towerinfosOH3->encode_key(channel);
0450     int this_etabin = towerinfosOH3->getTowerEtaBin(key);
0451     int this_phibin = towerinfosOH3->getTowerPhiBin(key);
0452     TowerInfo *tower = towerinfosOH3->get_tower_at_channel(channel);
0453     float this_E = tower->get_energy();
0454     int this_isBad = tower->get_isHot() || tower->get_isNoCalib() || tower->get_isNotInstr() || tower->get_isBadChi2();
0455     _OHCAL_ISBAD[this_etabin][this_phibin] = this_isBad;
0456     if (!this_isBad)
0457     { // just in case since all energy is summed
0458       _OHCAL_E[this_etabin][this_phibin] += this_E;
0459     }
0460
0461   }
0462
0463
0464   // first, calculate flow: Psi2 & v2, if enabled
0465
0466   //_Psi2 is left as 0
0467   // since _v2 is derived from _Psi2, initializing _Psi2 to NAN will set _v2 = NAN
0468   //_is_flow_failure tags events where _Psi2 & _v2 are set to 0 because there are no strips for flow
0469   // and when sEPD _Psi2 has no determined _Psi2 because the event is outside +/- z = 60cm
0470   _Psi2 = 0;
0471   _v2 = 0;
0472   _nStrips = 0;
0473   _is_flow_failure = false;
0474
0475
0476   if (_do_flow == 0)
0477   {
0478     if (Verbosity() > 0)
0479     {
0480       std::cout << "DetermineTowerBackground::process_event: flow not enabled, setting Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi ) , v2 = " << _v2 << std::endl;
0481     }
0482   }
0483
0484   if ( _do_flow >= 1 )
0485   {
0486
0487     if (Verbosity() > 0)
0488     {
0489       std::cout << "DetermineTowerBackground::process_event: flow enabled, calculating flow..." << std::endl;
0490     }
0491
0492     // phi weights are set to 1.0 by default
0493     _EMCAL_PHI_WEIGHTS.assign(_HCAL_NPHI, 1.0);
0494     _IHCAL_PHI_WEIGHTS.assign(_HCAL_NPHI, 1.0);
0495     _OHCAL_PHI_WEIGHTS.assign(_HCAL_NPHI, 1.0);
0496
0497     // copy the set of included eta strips to a new set for exclusion
0498     std::set<int> AVAILIBLE_ETA_STRIPS_CEMC = EtaStripsAvailbleForFlow;
0499     std::set<int> AVAILIBLE_ETA_STRIPS_IHCAL = EtaStripsAvailbleForFlow;
0500     std::set<int> AVAILIBLE_ETA_STRIPS_OHCAL = EtaStripsAvailbleForFlow;
0501
0502
0503     if ( _do_reweight )
0504     {
0505
0506       _reweight_failed = false; // flag to indicate if reweighting failed, i.e. all eta strips for a phi bin are excluded
0507       // rather than excluding full eta strips, we will reweight the phi bins according to the number of available eta strips
0508       // this is done by counting the number of eta strips which are still available for flow determination
0509       if (Verbosity() > 0)
0510       {
0511         std::cout << "DetermineTowerBackground::process_event: reweighting enabled, checking for bad towers in avialible eta strips..." << std::endl;
0512       }
0513
0514       // initialize the maximum number of eta bins per phi bin to be the total number of eta strips available (after removing the seeds)
0515       // loop over all phi bins
0516       for ( int phi = 0; phi < _HCAL_NPHI; phi++ )
0517       {
0518
0519         // initialize the maximum number of eta bins for this phi bin
0520         //  to be the total number of eta strips available (after removing the seeds)
0521         int EMCAL_MAX_TOWERS_THIS_PHI = MaxEtaBinsWithoutSeeds;
0522         int IHCAL_MAX_TOWERS_THIS_PHI = MaxEtaBinsWithoutSeeds;
0523         int OHCAL_MAX_TOWERS_THIS_PHI = MaxEtaBinsWithoutSeeds;
0524
0525         // loop over only the eta strips which are still available for flow determination
0526         // already exclude strips don't get checked since they are not in the set EtaStripsAvailbleForFlow
0527         for ( const auto &eta : EtaStripsAvailbleForFlow )
0528         {
0529
0530           EMCAL_MAX_TOWERS_THIS_PHI-= _EMCAL_ISBAD[eta][phi]; // decrement the possible count for this phi bin
0531           IHCAL_MAX_TOWERS_THIS_PHI-= _IHCAL_ISBAD[eta][phi]; // decrement the possible count for this phi bin
0532           OHCAL_MAX_TOWERS_THIS_PHI-= _OHCAL_ISBAD[eta][phi]; // decrement the possible count for this phi bin
0533           if ( Verbosity() > 10 )
0534           {
0535             if ( _EMCAL_ISBAD[eta][phi] )
0536             {
0537               std::cout << "DetermineTowerBackground::process_event: --> found bad tower in EMCAL at ieta / iphi = " << eta << " / " << phi << std::endl;
0538             }
0539             if ( _IHCAL_ISBAD[eta][phi] )
0540             {
0541               std::cout << "DetermineTowerBackground::process_event: --> found bad tower in IHCAL at ieta / iphi = " << eta << " / " << phi << std::endl;
0542             }
0543             if ( _OHCAL_ISBAD[eta][phi] )
0544             {
0545               std::cout << "DetermineTowerBackground::process_event: --> found bad tower in OHCAL at ieta / iphi = " << eta << " / " << phi << std::endl;
0546             }
0547           }
0548
0549         } // end loop over eta strips
0550
0551         if (Verbosity() > 1 )
0552         {
0553           std::cout << "DetermineTowerBackground::process_event: --> after checking for bad towers, EMCAL / IHCAL / OHCAL max eta strips for phi = "
0554             << phi << " are: " << EMCAL_MAX_TOWERS_THIS_PHI << " / " << IHCAL_MAX_TOWERS_THIS_PHI << " / " << OHCAL_MAX_TOWERS_THIS_PHI << std::endl;
0555         }
0556
0557         // update the phi weights for this phi bin
0558         if ( EMCAL_MAX_TOWERS_THIS_PHI > 0 )
0559         {
0560           _EMCAL_PHI_WEIGHTS[phi] = static_cast<float>(MaxEtaBinsWithoutSeeds) / static_cast<float>(EMCAL_MAX_TOWERS_THIS_PHI);
0561           if (Verbosity() > 0)
0562           {
0563             std::cout << "DetermineTowerBackground::process_event: --> setting EMCAL phi weight for phi = " << phi << " to " << _EMCAL_PHI_WEIGHTS[phi] << std::endl;
0564           }
0565         }
0566         else
0567         {
0568           // all the eta strips for this phi bin are excluded (this shouldn't happen)
0569           _EMCAL_PHI_WEIGHTS[phi] = 1.0;
0570           if (Verbosity() > 0)
0571           {
0572             std::cout << "DetermineTowerBackground::process_event: --> WARNING: all eta strips for EMCAL phi = " << phi << " are excluded, setting weight to 1.0" << std::endl;
0573             std::cout << "DeterminingTowerBackground::process_event: --> Defaulting to unweighted flow determination for this event." << std::endl;
0574           }
0575           _reweight_failed = true;
0576         }
0577         if ( IHCAL_MAX_TOWERS_THIS_PHI > 0 )
0578         {
0579           _IHCAL_PHI_WEIGHTS[phi] = static_cast<float>(MaxEtaBinsWithoutSeeds) / static_cast<float>(IHCAL_MAX_TOWERS_THIS_PHI);
0580           if (Verbosity() > 0)
0581           {
0582             std::cout << "DetermineTowerBackground::process_event: --> setting IHCAL phi weight for phi = " << phi << " to " << _IHCAL_PHI_WEIGHTS[phi] << std::endl;
0583           }
0584         }
0585         else
0586         {
0587           // all the eta strips for this phi bin are excluded (this shouldn't happen)
0588           _IHCAL_PHI_WEIGHTS[phi] = 1.0;
0589           if (Verbosity() > 0)
0590           {
0591             std::cout << "DetermineTowerBackground::process_event: --> WARNING: all eta strips for IHCAL phi = " << phi << " are excluded, setting weight to 1.0" << std::endl;
0592             std::cout << "DeterminingTowerBackground::process_event: --> Defaulting to unweighted flow determination for this event." << std::endl;
0593           }
0594           _reweight_failed = true;
0595
0596         }
0597         if ( OHCAL_MAX_TOWERS_THIS_PHI > 0 )
0598         {
0599           _OHCAL_PHI_WEIGHTS[phi] = static_cast<float>(MaxEtaBinsWithoutSeeds) / OHCAL_MAX_TOWERS_THIS_PHI;
0600           if (Verbosity() > 0)
0601           {
0602             std::cout << "DetermineTowerBackground::process_event: --> setting OHCAL phi weight for phi = " << phi << " to " << _OHCAL_PHI_WEIGHTS[phi] << std::endl;
0603           }
0604         }
0605         else
0606         {
0607           // all the eta strips for this phi bin are excluded (this shouldn't happen)
0608           _OHCAL_PHI_WEIGHTS[phi] = 1.0;
0609           if (Verbosity() > 0)
0610           {
0611             std::cout << "DetermineTowerBackground::process_event: --> WARNING: all eta strips for OHCAL phi = " << phi << " are excluded, setting weight to 1.0" << std::endl;
0612             std::cout << "DeterminingTowerBackground::process_event: --> Defaulting to unweighted flow determination for this event." << std::endl;
0613           }
0614           _reweight_failed = true;
0615
0616         }
0617       } // end loop over phi bins
0618
0619     }
0620
0621     if  (!_do_reweight || _reweight_failed )
0622     {
0623       // if reweighting is not enabled, we will exclude the eta strips which have bad towers in them
0624       if (Verbosity() > 0)
0625       {
0626         std::cout << "DetermineTowerBackground::process_event: reweighting not enabled, checking for bad towers in avialible eta strips..." << std::endl;
0627       }
0628       // loop over all available eta strips
0629       for ( const auto &eta : EtaStripsAvailbleForFlow )
0630       {
0631         if (Verbosity() > 2)
0632         {
0633           std::cout << "DetermineTowerBackground::process_event: checking for bad towers in eta strip " << eta << std::endl;
0634         }
0635         // get the number of bad phi towers within this eta strip
0636         // only look at the eta strips which are still available for flow determination which are in the
0637         // set EtaStripsAvailbleForFlow
0638         int bad_phis_int_this_eta_EMCAL = std::count(_EMCAL_ISBAD[eta].begin(), _EMCAL_ISBAD[eta].end(), 1); // count bad towers in this eta strip
0639         int bad_phis_int_this_eta_IHCAL = std::count(_IHCAL_ISBAD[eta].begin(), _IHCAL_ISBAD[eta].end(), 1); // count bad towers in this eta strip
0640         int bad_phis_int_this_eta_OHCAL = std::count(_OHCAL_ISBAD[eta].begin(), _OHCAL_ISBAD[eta].end(), 1); // count bad towers in this eta strip
0641         if (Verbosity() > 3)
0642         {
0643           std::cout << "DetermineTowerBackground::process_event: --> found " << bad_phis_int_this_eta_EMCAL << " bad towers in EMCAL, "
0644             << bad_phis_int_this_eta_IHCAL << " in IHCAL, and " << bad_phis_int_this_eta_OHCAL << " in OHCAL for eta strip " << eta << std::endl;
0645         }
0646         // we will exclude this eta strip if there are any bad towers in it
0647         if ( bad_phis_int_this_eta_EMCAL > 0 )
0648         {
0649           if (Verbosity() > 2)
0650           {
0651             std::cout << "DetermineTowerBackground::process_event: --> excluding EMCAL eta strip " << eta << " due to " << bad_phis_int_this_eta_EMCAL << " bad towers" << std::endl;
0652           }
0653           // remove this eta strip from the set of available eta strips
0654           AVAILIBLE_ETA_STRIPS_CEMC.erase(eta);
0655         }
0656         else
0657         {
0658           if (Verbosity() > 4)
0659           {
0660             std::cout << "DetermineTowerBackground::process_event: --> EMCAL eta strip " << eta << " has no excluded towers and can be used for flow determination " << std::endl;
0661           }
0662         }
0663
0664         if ( bad_phis_int_this_eta_IHCAL > 0 )
0665         {
0666           if (Verbosity() > 2)
0667           {
0668             std::cout << "DetermineTowerBackground::process_event: --> excluding IHCAL eta strip " << eta << " due to " << bad_phis_int_this_eta_IHCAL << " bad towers" << std::endl;
0669           }
0670           // remove this eta strip from the set of available eta strips
0671           AVAILIBLE_ETA_STRIPS_IHCAL.erase(eta);
0672         }
0673         else
0674         {
0675           if (Verbosity() > 4)
0676           {
0677             std::cout << "DetermineTowerBackground::process_event: --> IHCAL eta strip " << eta << " has no excluded towers and can be used for flow determination " << std::endl;
0678           }
0679         }
0680
0681         if ( bad_phis_int_this_eta_OHCAL > 0 )
0682         {
0683           if (Verbosity() > 2)
0684           {
0685             std::cout << "DetermineTowerBackground::process_event: --> excluding OHCAL eta strip " << eta << " due to " << bad_phis_int_this_eta_OHCAL << " bad towers" << std::endl;
0686           }
0687           // remove this eta strip from the set of available eta strips
0688           AVAILIBLE_ETA_STRIPS_OHCAL.erase(eta);
0689         }
0690         else
0691         {
0692           if (Verbosity() > 4)
0693           {
0694             std::cout << "DetermineTowerBackground::process_event: --> OHCAL eta strip " << eta << " has no excluded towers and can be used for flow determination " << std::endl;
0695           }
0696         }
0697
0698       } // end loop over eta strips
0699       if (Verbosity() > 0)
0700       {
0701         std::cout << "DetermineTowerBackground::process_event: after checking for bad towers, available EMCAL eta strips = " << AVAILIBLE_ETA_STRIPS_CEMC.size()
0702           << ", IHCAL eta strips = " << AVAILIBLE_ETA_STRIPS_IHCAL.size()
0703           << ", OHCAL eta strips = " << AVAILIBLE_ETA_STRIPS_OHCAL.size() << std::endl;
0704       }
0705     }
0706
0707     int nStripsAvailableForFlow = AVAILIBLE_ETA_STRIPS_CEMC.size() + AVAILIBLE_ETA_STRIPS_IHCAL.size() + AVAILIBLE_ETA_STRIPS_OHCAL.size();
0708     int nStripsUnavailableForFlow = (_HCAL_NETA*3) - nStripsAvailableForFlow;
0709     if (Verbosity() > 0)
0710     {
0711       std::cout << "DetermineTowerBackground::process_event: # of strips (summed over layers) available / unavailable for flow determination: " << nStripsAvailableForFlow << " / " << nStripsUnavailableForFlow << std::endl;
0712     }
0713
0714     if ( nStripsAvailableForFlow > 0 )
0715     {
0716
0717
0718       _nStrips = nStripsAvailableForFlow;
0719
0720       // update the full calorimeter flow vectors
0721       _FULLCALOFLOW_PHI_E.assign(_HCAL_NPHI, 0.0);
0722       _FULLCALOFLOW_PHI_VAL.assign(_HCAL_NPHI, 0.0);
0723
0724       // flow determination
0725       float Q_x = 0;
0726       float Q_y = 0;
0727       float sum_E = 0;
0728       for (int phi = 0; phi < _HCAL_NPHI; phi++)
0729       {
0730         _FULLCALOFLOW_PHI_VAL[phi] = geomIH->get_phicenter(phi);
0731         // loop over the available eta strips for each layer
0732
0733         for (const auto &eta : AVAILIBLE_ETA_STRIPS_CEMC)
0734         {
0735           _FULLCALOFLOW_PHI_E[phi] += _EMCAL_E[eta][phi] * _EMCAL_PHI_WEIGHTS[phi]; // if reweighting is enabled, the weights are applied, if not, they are 1.0
0736         }
0737         for (const auto &eta : AVAILIBLE_ETA_STRIPS_IHCAL)
0738         {
0739           _FULLCALOFLOW_PHI_E[phi] += _IHCAL_E[eta][phi] * _IHCAL_PHI_WEIGHTS[phi]; // if reweighting is enabled, the weights are applied, if not, they are 1.0
0740         }
0741         for (const auto &eta : AVAILIBLE_ETA_STRIPS_OHCAL)
0742         {
0743           _FULLCALOFLOW_PHI_E[phi] += _OHCAL_E[eta][phi] * _OHCAL_PHI_WEIGHTS[phi]; // if reweighting is enabled, the weights are applied, if not, they are 1.0
0744         }
0745
0746         // sum up the energy in this phi bin
0747         Q_x += _FULLCALOFLOW_PHI_E[phi] * cos(2 * _FULLCALOFLOW_PHI_VAL[phi]);
0748         Q_y += _FULLCALOFLOW_PHI_E[phi] * sin(2 * _FULLCALOFLOW_PHI_VAL[phi]);
0749         sum_E += _FULLCALOFLOW_PHI_E[phi];
0750
0751       }
0752
0753       if (_do_flow == 1)
0754       { // Calo event plane
0755         _Psi2 = std::atan2(Q_y, Q_x) / 2.0;
0756       }
0757       else if (_do_flow == 2)
0758       { // HIJING truth flow extraction
0759         PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0760
0761         if (!truthinfo)
0762         {
0763           std::cout << "DetermineTowerBackground::process_event: FATAL , G4TruthInfo does not exist , cannot extract truth flow with do_flow = " << _do_flow << std::endl;
0764           return -1;
0765         }
0766
0767         PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
0768
0769         float Hijing_Qx = 0;
0770         float Hijing_Qy = 0;
0771
0772         for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
0773         {
0774           PHG4Particle *g4particle = iter->second;
0775
0776           if (truthinfo->isEmbeded(g4particle->get_track_id()) != 0)
0777           {
0778             continue;
0779           }
0780
0781           TLorentzVector t;
0782           t.SetPxPyPzE(g4particle->get_px(), g4particle->get_py(), g4particle->get_pz(), g4particle->get_e());
0783
0784           float truth_pt = t.Pt();
0785           if (truth_pt < 0.4)
0786           {
0787             continue;
0788           }
0789           float truth_eta = t.Eta();
0790           if (std::fabs(truth_eta) > 1.1)
0791           {
0792             continue;
0793           }
0794           float truth_phi = t.Phi();
0795           int truth_pid = g4particle->get_pid();
0796
0797           if (Verbosity() > 10)
0798           {
0799             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;
0800           }
0801
0802           Hijing_Qx += truth_pt * std::cos(2 * truth_phi);
0803           Hijing_Qy += truth_pt * std::sin(2 * truth_phi);
0804         }
0805
0806         _Psi2 = std::atan2(Hijing_Qy, Hijing_Qx) / 2.0;
0807
0808         if (Verbosity() > 0)
0809         {
0810           std::cout << "DetermineTowerBackground::process_event: flow extracted from Hijing truth particles, setting Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi ) " << std::endl;
0811         }
0812       }
0813       else if (_do_flow == 3)
0814       { // sEPD event plane extraction
0815         // get event plane map
0816         EventplaneinfoMap *epmap = findNode::getClass<EventplaneinfoMap>(topNode, "EventplaneinfoMap");
0817         if (!epmap)
0818         {
0819           std::cout << "DetermineTowerBackground::process_event: FATAL, EventplaneinfoMap does not exist, cannot extract sEPD flow with do_flow = " << _do_flow << std::endl;
0820           exit(-1);
0821         }
0822         if (!(epmap->empty()))
0823         {
0824           auto *EPDNS = epmap->get(EventplaneinfoMap::sEPDNS);
0825           _Psi2 = EPDNS->get_shifted_psi(2);
0826         }
0827         else
0828         {
0829           _is_flow_failure = true;
0830           _Psi2 = 0;
0831         }
0832
0833         if (Verbosity() > 0)
0834         {
0835           std::cout << "DetermineTowerBackground::process_event: flow extracted from sEPD, setting Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi ) " << std::endl;
0836         }
0837
0838       }
0839
0840       if (std::isnan(_Psi2) || std::isinf(_Psi2))
0841       {
0842         _Psi2 = 0;
0843         _is_flow_failure = true;
0844         if (Verbosity() > 0)
0845         {
0846           std::cout << "DetermineTowerBackground::process_event: sEPD event plane extraction failed, setting Psi2 = 0" << std::endl;
0847         }
0848       }
0849
0850
0851       _v2 = 0;
0852       for (int phi = 0; phi < _HCAL_NPHI; phi++)
0853       {
0854         _v2 += ( _FULLCALOFLOW_PHI_E[phi] * std::cos(2 * (_FULLCALOFLOW_PHI_VAL[phi] - _Psi2)) );
0855       }
0856
0857
0858       // avoid nans in v2
0859       if (sum_E > 0)
0860       {
0861         _v2 /= sum_E;
0862       }
0863       else
0864       {
0865         _v2 = 0;
0866       }
0867
0868       if (Verbosity() > 0)
0869       {
0870         std::cout << "DetermineTowerBackground::process_event: unnormalized Q vector (Qx, Qy) = ( " << Q_x << ", " << Q_y << " ) with Sum E_i = " << sum_E << std::endl;
0871         std::cout << "DetermineTowerBackground::process_event: Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi " << (_do_flow == 2 ? "from Hijing " : "") << ") , v2 = " << _v2 << " ( using " << _nStrips << " ) " << std::endl;
0872       }
0873     }
0874     else
0875     {
0876       _Psi2 = 0;
0877       _v2 = 0;
0878       _nStrips = 0;
0879       _is_flow_failure = true;
0880       if (Verbosity() > 0)
0881       {
0882         std::cout << "DetermineTowerBackground::process_event: no full strips available for flow modulation, setting v2 and Psi = 0" << std::endl;
0883       }
0884
0885     }
0886
0887
0888     if (Verbosity() > 0)
0889     {
0890       std::cout << "DetermineTowerBackground::process_event: flow extraction successful, Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi ) , v2 = " << _v2 << std::endl;
0891     }
0892   }  // if do flow
0893
0894
0895   // now calculate energy densities...
0896   _nTowers = 0;  // store how many towers were used to determine bkg
0897
0898   // starting with the EMCal first...
0899   for (int layer = 0; layer < 3; layer++)
0900   {
0901     int local_max_eta = _HCAL_NETA;
0902     int local_max_phi = _HCAL_NPHI;
0903
0904     for (int eta = 0; eta < local_max_eta; eta++)
0905     {
0906       float total_E = 0;
0907       int total_tower = 0;
0908
0909       for (int phi = 0; phi < local_max_phi; phi++)
0910       {
0911         float this_eta = geomIH->get_etacenter(eta);
0912         float this_phi = geomIH->get_phicenter(phi);
0913
0914         bool isExcluded = false;
0915
0916         // nobody can understand these nested ?s, which just makes this error prone and a maintenance headache
0917         //        float my_E = (layer == 0 ? _EMCAL_E[eta][phi] : (layer == 1 ? _IHCAL_E[eta][phi] : _OHCAL_E[eta][phi]));
0918         //        int this_isBad = (layer == 0 ? _EMCAL_ISBAD[eta][phi] : (layer == 1 ? _IHCAL_ISBAD[eta][phi] : _OHCAL_ISBAD[eta][phi]));
0919         //  please use the following if/else if/else
0920         float my_E;
0921         int this_isBad;
0922         if (layer == 0)
0923         {
0924           my_E = _EMCAL_E[eta][phi];
0925           this_isBad = _EMCAL_ISBAD[eta][phi];
0926         }
0927         else if (layer == 1)
0928         {
0929           my_E = _IHCAL_E[eta][phi];
0930           this_isBad = _IHCAL_ISBAD[eta][phi];
0931         }
0932         else
0933         {
0934           my_E = _OHCAL_E[eta][phi];
0935           this_isBad = _OHCAL_ISBAD[eta][phi];
0936         }
0937         // if the tower is masked (energy identically zero), exclude it
0938         if (this_isBad)
0939         {
0940           isExcluded = true;
0941           if (Verbosity() > 10)
0942           {
0943             std::cout << " tower in layer " << layer << " at eta / phi = " << this_eta << " / " << this_phi << " with E = " << my_E << " excluded due to masking" << std::endl;
0944           }
0945         }
0946         for (unsigned int iseed = 0; iseed < _seed_eta.size(); iseed++)
0947         {
0948           float deta = this_eta - _seed_eta[iseed];
0949           float dphi = this_phi - _seed_phi[iseed];
0950           if (dphi > M_PI)
0951           {
0952             dphi -= 2 * M_PI;
0953           }
0954           if (dphi < -M_PI)
0955           {
0956             dphi += 2 * M_PI;
0957           }
0958           float dR = sqrt(pow(deta, 2) + pow(dphi, 2));
0959           if (dR < 0.4)
0960           {
0961             isExcluded = true;
0962             if (Verbosity() > 10)
0963             {
0964               std::cout << " setting excluded mark from seed at eta / phi = " << _seed_eta[iseed] << " / " << _seed_phi[iseed] << std::endl;
0965             }
0966           }
0967         }
0968
0969         float modulation_factor =  1 + 2 * _v2 * std::cos(2 * (this_phi - _Psi2));
0970         if (!isExcluded && modulation_factor > 0)
0971         {
0972           if (layer == 0)
0973           {
0974             total_E += _EMCAL_E[eta][phi] / modulation_factor;
0975           }
0976           if (layer == 1)
0977           {
0978             total_E += _IHCAL_E[eta][phi] / modulation_factor;
0979           }
0980           if (layer == 2)
0981           {
0982             total_E += _OHCAL_E[eta][phi] / modulation_factor;
0983           }
0984           total_tower++;  // towers in this eta range & layer
0985           _nTowers++;     // towers in entire calorimeter
0986         }
0987         else
0988         {
0989           if (Verbosity() > 10)
0990           {
0991             std::cout << " tower at eta / phi = " << this_eta << " / " << this_phi << " with E = " << total_E << " excluded due to seed " << std::endl;
0992           }
0993         }
0994       }
0995
0996       std::pair<float, float> etabounds = geomIH->get_etabounds(eta);
0997       std::pair<float, float> phibounds = geomIH->get_phibounds(0);
0998
0999       float deta = etabounds.second - etabounds.first;
1000       float dphi = phibounds.second - phibounds.first;
1001       float total_area = total_tower * deta * dphi;
1002
1003       if ( total_tower > 0 )
1004       {
1005         _UE[layer].at(eta) = total_E / total_tower; // calculate the UE density
1006       }
1007       else
1008       {
1009         if (Verbosity() > 0)
1010         {
1011           std::cout << "DetermineTowerBackground::process_event: WARNING, no towers in layer " << layer << " / eta " << eta << ", setting UE density to 0" << std::endl;
1012         }
1013
1014         _UE[layer].at(eta) = 0; // no towers, no UE
1015       }
1016
1017       if (Verbosity() > 3)
1018       {
1019         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;
1020       }
1021     }
1022   }
1023
1024   if (Verbosity() > 0)
1025   {
1026     for (int layer = 0; layer < 3; layer++)
1027     {
1028       std::cout << "DetermineTowerBackground::process_event: summary UE in layer " << layer << " : ";
1029       for (int eta = 0; eta < _HCAL_NETA; eta++)
1030       {
1031         std::cout << _UE[layer].at(eta) << " , ";
1032       }
1033       std::cout << std::endl;
1034     }
1035   }
1036
1037   FillNode(topNode);
1038
1039   if (Verbosity() > 0)
1040   {
1041     std::cout << "DetermineTowerBackground::process_event: exiting" << std::endl;
1042   }
1043
1044   return Fun4AllReturnCodes::EVENT_OK;
1045 }
1046
1047 int DetermineTowerBackground::CreateNode(PHCompositeNode *topNode)
1048 {
1049   PHNodeIterator iter(topNode);
1050
1051   // Looking for the DST node
1052   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
1053   if (!dstNode)
1054   {
1055     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
1056     return Fun4AllReturnCodes::ABORTRUN;
1057   }
1058
1059   // store the jet background stuff under a sub-node directory
1060   PHCompositeNode *bkgNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "JETBACKGROUND"));
1061   if (!bkgNode)
1062   {
1063     bkgNode = new PHCompositeNode("JETBACKGROUND");
1064     dstNode->addNode(bkgNode);
1065   }
1066
1067   // create the TowerBackground node...
1068   TowerBackground *towerbackground = findNode::getClass<TowerBackground>(topNode, _backgroundName);
1069   if (!towerbackground)
1070   {
1071     towerbackground = new TowerBackgroundv1();
1072     PHIODataNode<PHObject> *bkgDataNode = new PHIODataNode<PHObject>(towerbackground, _backgroundName, "PHObject");
1073     bkgNode->addNode(bkgDataNode);
1074   }
1075   else
1076   {
1077     std::cout << PHWHERE << "::ERROR - " << _backgroundName << " pre-exists, but should not" << std::endl;
1078     // exit(-1);
1079   }
1080
1081   return Fun4AllReturnCodes::EVENT_OK;
1082 }
1083
1084 void DetermineTowerBackground::FillNode(PHCompositeNode *topNode)
1085 {
1086   TowerBackground *towerbackground = findNode::getClass<TowerBackground>(topNode, _backgroundName);
1087   if (!towerbackground)
1088   {
1089     std::cout << " ERROR -- can't find TowerBackground node after it should have been created" << std::endl;
1090     return;
1091   }
1092
1093   towerbackground->set_UE(0, _UE[0]);
1094   towerbackground->set_UE(1, _UE[1]);
1095   towerbackground->set_UE(2, _UE[2]);
1096
1097   towerbackground->set_v2(_v2);
1098
1099   towerbackground->set_Psi2(_Psi2);
1100
1101   towerbackground->set_nStripsUsedForFlow(_nStrips);
1102
1103   towerbackground->set_nTowersUsedForBkg(_nTowers);
1104
1105   towerbackground->set_flow_failure_flag(_is_flow_failure);
1106
1107   return;
1108 }
1109
1110
1111
1112
1113