Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:17:43 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!