Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-12-16 09:20:03

0001 #include "SubtractTowers.h"
0002 
0003 #include "TowerBackground.h"
0004 
0005 // sPHENIX includes
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/RawTowerv1.h>
0012 
0013 #include <calobase/TowerInfo.h>
0014 #include <calobase/TowerInfoContainer.h>
0015 
0016 #include <fun4all/Fun4AllReturnCodes.h>
0017 #include <fun4all/SubsysReco.h>
0018 
0019 #include <phool/PHCompositeNode.h>
0020 #include <phool/PHIODataNode.h>
0021 #include <phool/PHNode.h>
0022 #include <phool/PHNodeIterator.h>
0023 #include <phool/PHObject.h>
0024 #include <phool/getClass.h>
0025 #include <phool/phool.h>
0026 
0027 // standard includes
0028 #include <cmath>
0029 #include <cstdlib>
0030 #include <iostream>
0031 #include <map>
0032 #include <utility>
0033 #include <vector>
0034 
0035 SubtractTowers::SubtractTowers(const std::string &name)
0036   : SubsysReco(name)
0037 {
0038 }
0039 
0040 int SubtractTowers::InitRun(PHCompositeNode *topNode)
0041 {
0042   CreateNode(topNode);
0043 
0044   return Fun4AllReturnCodes::EVENT_OK;
0045 }
0046 
0047 int SubtractTowers::process_event(PHCompositeNode *topNode)
0048 {
0049   if (Verbosity() > 0)
0050   {
0051     std::cout << "SubtractTowers::process_event: entering, with _use_flow_modulation = " << _use_flow_modulation << std::endl;
0052   }
0053 
0054   // pull out the tower containers and geometry objects at the start
0055   RawTowerContainer *towersEM3 = nullptr;
0056   RawTowerContainer *towersIH3 = nullptr;
0057   RawTowerContainer *towersOH3 = nullptr;
0058   TowerInfoContainer *towerinfosEM3 = nullptr;
0059   TowerInfoContainer *towerinfosIH3 = nullptr;
0060   TowerInfoContainer *towerinfosOH3 = nullptr;
0061 
0062   if (m_use_towerinfo)
0063   {
0064     EMTowerName = m_towerNodePrefix + "_CEMC_RETOWER";
0065     IHTowerName = m_towerNodePrefix + "_HCALIN";
0066     OHTowerName = m_towerNodePrefix + "_HCALOUT";
0067     towerinfosEM3 = findNode::getClass<TowerInfoContainer>(topNode, EMTowerName);
0068     towerinfosIH3 = findNode::getClass<TowerInfoContainer>(topNode, IHTowerName);
0069     towerinfosOH3 = findNode::getClass<TowerInfoContainer>(topNode, OHTowerName);
0070 
0071     if (Verbosity() > 0)
0072     {
0073       std::cout << "SubtractTowers::process_event: " << towerinfosEM3->size() << EMTowerName << " towers" << std::endl;
0074       std::cout << "SubtractTowers::process_event: " << towerinfosIH3->size() << IHTowerName << " towers" << std::endl;
0075       std::cout << "SubtractTowers::process_event: " << towerinfosOH3->size() << OHTowerName << " towers" << std::endl;
0076     }
0077   }
0078   else
0079   {
0080     towersEM3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER");
0081     towersIH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
0082     towersOH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
0083 
0084     if (Verbosity() > 0)
0085     {
0086       std::cout << "SubtractTowers::process_event: " << towersEM3->size() << " TOWER_CALIB_CEMC_RETOWER towers" << std::endl;
0087       std::cout << "SubtractTowers::process_event: " << towersIH3->size() << " TOWER_CALIB_HCALIN towers" << std::endl;
0088       std::cout << "SubtractTowers::process_event: " << towersOH3->size() << " TOWER_CALIB_HCALOUT towers" << std::endl;
0089     }
0090   }
0091 
0092   RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0093   RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0094 
0095   // these should have already been created during InitRun()
0096   RawTowerContainer *emcal_towers = nullptr;
0097   RawTowerContainer *ihcal_towers = nullptr;
0098   RawTowerContainer *ohcal_towers = nullptr;
0099   TowerInfoContainer *emcal_towerinfos = nullptr;
0100   TowerInfoContainer *ihcal_towerinfos = nullptr;
0101   TowerInfoContainer *ohcal_towerinfos = nullptr;
0102   if (m_use_towerinfo)
0103   {
0104     EMTowerName = m_towerNodePrefix + "_CEMC_RETOWER_SUB1";
0105     IHTowerName = m_towerNodePrefix + "_HCALIN_SUB1";
0106     OHTowerName = m_towerNodePrefix + "_HCALOUT_SUB1";
0107     emcal_towerinfos = findNode::getClass<TowerInfoContainer>(topNode, EMTowerName);
0108     ihcal_towerinfos = findNode::getClass<TowerInfoContainer>(topNode, IHTowerName);
0109     ohcal_towerinfos = findNode::getClass<TowerInfoContainer>(topNode, OHTowerName);
0110   
0111     if (Verbosity() > 0)
0112       {
0113     std::cout << "SubtractTowers::process_event: starting with " << emcal_towerinfos->size() << EMTowerName << " towers" << std::endl;
0114     std::cout << "SubtractTowers::process_event: starting with " << ihcal_towerinfos->size() << IHTowerName << " towers" << std::endl;
0115     std::cout << "SubtractTowers::process_event: starting with " << ohcal_towerinfos->size() << OHTowerName << " towers" << std::endl;
0116       }
0117   }
0118 
0119   else
0120   {
0121     emcal_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER_SUB1");
0122     ihcal_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN_SUB1");
0123     ohcal_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT_SUB1");
0124     if (Verbosity() > 0)
0125     {
0126       std::cout << "SubtractTowers::process_event: starting with " << emcal_towers->size() << " TOWER_CALIB_CEMC_RETOWER_SUB1 towers" << std::endl;
0127       std::cout << "SubtractTowers::process_event: starting with " << ihcal_towers->size() << " TOWER_CALIB_HCALIN_SUB1 towers" << std::endl;
0128       std::cout << "SubtractTowers::process_event: starting with " << ohcal_towers->size() << " TOWER_CALIB_HCALOUT_SUB1 towers" << std::endl;
0129     }
0130   }
0131 
0132   TowerBackground *towerbackground;
0133   if (m_use_towerinfo)
0134   {
0135     towerbackground = findNode::getClass<TowerBackground>(topNode, "TowerInfoBackground_Sub2");
0136   }
0137   else
0138   {
0139     towerbackground = findNode::getClass<TowerBackground>(topNode, "TowerBackground_Sub2");
0140   }
0141   // read these in to use, even if we don't use flow modulation in the subtraction
0142   float background_v2 = towerbackground->get_v2();
0143   float background_Psi2 = towerbackground->get_Psi2();
0144 
0145   // EMCal
0146 
0147   // replicate existing towers
0148   if (m_use_towerinfo)
0149   {
0150     unsigned int nchannels_em = towerinfosEM3->size();
0151     for (unsigned int channel = 0; channel < nchannels_em; channel++)
0152     {
0153       TowerInfo *tower = towerinfosEM3->get_tower_at_channel(channel);
0154       unsigned int towerkey = towerinfosEM3->encode_key(channel);
0155       int ieta = towerinfosEM3->getTowerEtaBin(towerkey);
0156       int iphi = towerinfosEM3->getTowerPhiBin(towerkey);
0157       float raw_energy = tower->get_energy();
0158       float UE = towerbackground->get_UE(0).at(ieta);
0159       if (_use_flow_modulation)
0160       {
0161         const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0162         float tower_phi = geomIH->get_tower_geometry(key)->get_phi();
0163         float modulation_factor = 1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2));
0164         modulation_factor = std::max(0.F, modulation_factor);
0165         UE = UE * modulation_factor;
0166       }
0167       float new_energy = raw_energy - UE;
0168       // if a tower is masked, leave it at zero
0169       if (tower->get_isHot() || tower->get_isNoCalib() || tower->get_isNotInstr() || tower->get_isBadChi2())
0170       {
0171         new_energy = 0;
0172       }
0173 
0174       emcal_towerinfos->get_tower_at_channel(channel)->set_time(tower->get_time());
0175       emcal_towerinfos->get_tower_at_channel(channel)->set_energy(new_energy);
0176 
0177       if (Verbosity() > 5)
0178       {
0179         std::cout << " SubtractTowers::process_event : EMCal tower at ieta / iphi = " << ieta << " / " << iphi << ", pre-sub / after-sub E = " << raw_energy << " / " << new_energy << std::endl;
0180       }
0181     }
0182   }
0183   else
0184   {
0185     RawTowerContainer::ConstRange begin_end_EM = towersEM3->getTowers();
0186     for (RawTowerContainer::ConstIterator rtiter = begin_end_EM.first; rtiter != begin_end_EM.second; ++rtiter)
0187     {
0188       RawTower *tower = rtiter->second;
0189 
0190       int this_etabin = tower->get_bineta();
0191       int this_phibin = tower->get_binphi();
0192       float this_E = tower->get_energy();
0193 
0194       RawTower *new_tower = new RawTowerv1();
0195       new_tower->set_energy(this_E);
0196       emcal_towers->AddTower(this_etabin, this_phibin, new_tower);
0197     }
0198 
0199     // now fill in additional towers with zero energy to fill out the full grid
0200     // but note: after retowering, all of these should already exist...
0201     for (int eta = 0; eta < geomIH->get_etabins(); eta++)
0202     {
0203       for (int phi = 0; phi < geomIH->get_phibins(); phi++)
0204       {
0205         if (!emcal_towers->getTower(eta, phi))
0206         {
0207           RawTower *new_tower = new RawTowerv1();
0208           new_tower->set_energy(0);
0209           emcal_towers->AddTower(eta, phi, new_tower);
0210         }
0211       }
0212     }
0213     // update towers for background subtraction...
0214     for (RawTowerContainer::ConstIterator rtiter = emcal_towers->getTowers().first; rtiter != emcal_towers->getTowers().second; ++rtiter)
0215     {
0216       RawTower *tower = rtiter->second;
0217       float raw_energy = tower->get_energy();
0218       float UE = towerbackground->get_UE(0).at(tower->get_bineta());
0219       if (_use_flow_modulation)
0220       {
0221         float tower_phi = geomIH->get_tower_geometry(tower->get_key())->get_phi();
0222         float modulation_factor = 1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2));
0223         modulation_factor = std::max(0.F, modulation_factor);
0224         UE = UE * modulation_factor;
0225       }
0226       float new_energy = raw_energy - UE;
0227       if (raw_energy == 0)
0228       {
0229         new_energy = 0;
0230       }
0231       tower->set_energy(new_energy);
0232       if (Verbosity() > 5)
0233       {
0234         std::cout << "SubtractTowers::process_event : EMCal tower at eta / phi = " << tower->get_bineta() << " / " << tower->get_binphi() << ", pre-sub / after-sub E = " << raw_energy << " / " << tower->get_energy() << std::endl;
0235       }
0236     }
0237   }
0238   // IHCal
0239   // replicate existing towers
0240   if (m_use_towerinfo)
0241   {
0242     unsigned int nchannels_ih = towerinfosIH3->size();
0243     for (unsigned int channel = 0; channel < nchannels_ih; channel++)
0244     {
0245       TowerInfo *tower = towerinfosIH3->get_tower_at_channel(channel);
0246       unsigned int towerkey = towerinfosIH3->encode_key(channel);
0247       int ieta = towerinfosIH3->getTowerEtaBin(towerkey);
0248       int iphi = towerinfosIH3->getTowerPhiBin(towerkey);
0249 
0250       float raw_energy = tower->get_energy();
0251       float UE = towerbackground->get_UE(1).at(ieta);
0252       if (_use_flow_modulation)
0253       {
0254         const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0255         float tower_phi = geomIH->get_tower_geometry(key)->get_phi();
0256         float modulation_factor = 1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2));
0257         modulation_factor = std::max(0.F, modulation_factor);
0258         UE = UE * modulation_factor;
0259       }
0260       float new_energy = raw_energy - UE;
0261       // if a tower is masked, leave it at zero
0262       if (tower->get_isHot() || tower->get_isNoCalib() || tower->get_isNotInstr() || tower->get_isBadChi2())
0263       {
0264         new_energy = 0;
0265       }
0266 
0267       ihcal_towerinfos->get_tower_at_channel(channel)->set_time(tower->get_time());
0268       ihcal_towerinfos->get_tower_at_channel(channel)->set_energy(new_energy);
0269       if (Verbosity() > 5)
0270       {
0271         std::cout << "SubtractTowers::process_event : IHCal tower at ieta / iphi = " << ieta << " / " << iphi << ", pre-sub / after-sub E = " << raw_energy << " / " << new_energy << std::endl;
0272       }
0273     }
0274   }
0275   else
0276   {
0277     RawTowerContainer::ConstRange begin_end_IH = towersIH3->getTowers();
0278     for (RawTowerContainer::ConstIterator rtiter = begin_end_IH.first; rtiter != begin_end_IH.second; ++rtiter)
0279     {
0280       RawTower *tower = rtiter->second;
0281       RawTowerGeom *tower_geom = geomIH->get_tower_geometry(tower->get_key());
0282 
0283       int this_etabin = geomIH->get_etabin(tower_geom->get_eta());
0284       int this_phibin = geomIH->get_phibin(tower_geom->get_phi());
0285       float this_E = tower->get_energy();
0286 
0287       RawTower *new_tower = new RawTowerv1();
0288       new_tower->set_energy(this_E);
0289       ihcal_towers->AddTower(this_etabin, this_phibin, new_tower);
0290     }
0291 
0292     // now fill in additional towers with zero energy to fill out the full grid
0293     for (int eta = 0; eta < geomIH->get_etabins(); eta++)
0294     {
0295       for (int phi = 0; phi < geomIH->get_phibins(); phi++)
0296       {
0297         if (!ihcal_towers->getTower(eta, phi))
0298         {
0299           RawTower *new_tower = new RawTowerv1();
0300           new_tower->set_energy(0);
0301           ihcal_towers->AddTower(eta, phi, new_tower);
0302         }
0303       }
0304     }
0305 
0306     // update towers for background subtraction...
0307     for (RawTowerContainer::ConstIterator rtiter = ihcal_towers->getTowers().first; rtiter != ihcal_towers->getTowers().second; ++rtiter)
0308     {
0309       RawTower *tower = rtiter->second;
0310       float raw_energy = tower->get_energy();
0311       float UE = towerbackground->get_UE(1).at(tower->get_bineta());
0312       if (_use_flow_modulation)
0313       {
0314         float tower_phi = geomIH->get_tower_geometry(tower->get_key())->get_phi();
0315         float modulation_factor = 1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2));
0316         modulation_factor = std::max(0.F, modulation_factor);
0317         UE = UE * modulation_factor;
0318       }
0319       float new_energy = raw_energy - UE;
0320       if (raw_energy == 0)
0321       {
0322         new_energy = 0;
0323       }
0324       tower->set_energy(new_energy);
0325     }
0326   }
0327   // OHCal
0328 
0329   // replicate existing towers
0330   if (m_use_towerinfo)
0331   {
0332     unsigned int nchannels_oh = towerinfosOH3->size();
0333     for (unsigned int channel = 0; channel < nchannels_oh; channel++)
0334     {
0335       TowerInfo *tower = towerinfosOH3->get_tower_at_channel(channel);
0336       unsigned int towerkey = towerinfosOH3->encode_key(channel);
0337       int ieta = towerinfosOH3->getTowerEtaBin(towerkey);
0338       int iphi = towerinfosOH3->getTowerPhiBin(towerkey);
0339       float raw_energy = tower->get_energy();
0340       float UE = towerbackground->get_UE(2).at(ieta);
0341       if (_use_flow_modulation)
0342       {
0343         const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0344         float tower_phi = geomOH->get_tower_geometry(key)->get_phi();
0345         float modulation_factor = 1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2));
0346         modulation_factor = std::max(0.F, modulation_factor);
0347         UE = UE * modulation_factor;
0348       }
0349       float new_energy = raw_energy - UE;
0350       // if a tower is masked, leave it at zero
0351       if (tower->get_isHot() || tower->get_isNoCalib() || tower->get_isNotInstr() || tower->get_isBadChi2())
0352       {
0353         new_energy = 0;
0354       }
0355 
0356       ohcal_towerinfos->get_tower_at_channel(channel)->set_time(tower->get_time());
0357       ohcal_towerinfos->get_tower_at_channel(channel)->set_energy(new_energy);
0358       if (Verbosity() > 5)
0359       {
0360         std::cout << "SubtractTowers::process_event : OHCal tower at ieta / iphi = " << ieta << " / " << iphi << ", pre-sub / after-sub E = " << raw_energy << " / " << new_energy << std::endl;
0361       }
0362     }
0363   }
0364   else
0365   {
0366     RawTowerContainer::ConstRange begin_end_OH = towersOH3->getTowers();
0367     for (RawTowerContainer::ConstIterator rtiter = begin_end_OH.first; rtiter != begin_end_OH.second; ++rtiter)
0368     {
0369       RawTower *tower = rtiter->second;
0370       RawTowerGeom *tower_geom = geomOH->get_tower_geometry(tower->get_key());
0371 
0372       int this_etabin = geomOH->get_etabin(tower_geom->get_eta());
0373       int this_phibin = geomOH->get_phibin(tower_geom->get_phi());
0374       float this_E = tower->get_energy();
0375 
0376       RawTower *new_tower = new RawTowerv1();
0377       new_tower->set_energy(this_E);
0378       ohcal_towers->AddTower(this_etabin, this_phibin, new_tower);
0379     }
0380 
0381     // now fill in additional towers with zero energy to fill out the full grid
0382     for (int eta = 0; eta < geomOH->get_etabins(); eta++)
0383     {
0384       for (int phi = 0; phi < geomOH->get_phibins(); phi++)
0385       {
0386         if (!ohcal_towers->getTower(eta, phi))
0387         {
0388           RawTower *new_tower = new RawTowerv1();
0389           new_tower->set_energy(0);
0390           ohcal_towers->AddTower(eta, phi, new_tower);
0391         }
0392       }
0393     }
0394 
0395     // update towers for background subtraction...
0396     for (RawTowerContainer::ConstIterator rtiter = ohcal_towers->getTowers().first; rtiter != ohcal_towers->getTowers().second; ++rtiter)
0397     {
0398       RawTower *tower = rtiter->second;
0399       float raw_energy = tower->get_energy();
0400       float UE = towerbackground->get_UE(2).at(tower->get_bineta());
0401       if (_use_flow_modulation)
0402       {
0403         float tower_phi = geomOH->get_tower_geometry(tower->get_key())->get_phi();
0404         float modulation_factor = 1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2));
0405         modulation_factor = std::max(0.F, modulation_factor);
0406         UE = UE * modulation_factor;
0407       }
0408       float new_energy = raw_energy - UE;
0409       if (raw_energy == 0)
0410       {
0411         new_energy = 0;
0412       }
0413       tower->set_energy(new_energy);
0414     }
0415   }
0416 
0417   if (Verbosity() > 0)
0418   {
0419     if (!m_use_towerinfo)
0420     {
0421       std::cout << "SubtractTowers::process_event: ending with " << emcal_towers->size() << " TOWER_CALIB_CEMC_RETOWER_SUB1 towers" << std::endl;
0422       std::cout << "SubtractTowers::process_event: ending with " << ihcal_towers->size() << " TOWER_CALIB_HCALIN_SUB1 towers" << std::endl;
0423       std::cout << "SubtractTowers::process_event: ending with " << ohcal_towers->size() << " TOWER_CALIB_HCALOUT_SUB1 towers" << std::endl;
0424     }
0425     else
0426     {
0427       std::cout << "SubtractTowers::process_event: ending with " << emcal_towerinfos->size() << m_towerNodePrefix << "_CEMC_RETOWER_SUB1 towers" << std::endl;
0428       std::cout << "SubtractTowers::process_event: ending with " << ihcal_towerinfos->size() << m_towerNodePrefix << "_HCALIN_SUB1 towers" << std::endl;
0429       std::cout << "SubtractTowers::process_event: ending with " << ohcal_towerinfos->size() << m_towerNodePrefix << "_HCALOUT_SUB1 towers" << std::endl;
0430     }
0431   }
0432 
0433   if (Verbosity() > 0)
0434   {
0435     std::cout << "SubtractTowers::process_event: exiting" << std::endl;
0436   }
0437 
0438   return Fun4AllReturnCodes::EVENT_OK;
0439 }
0440 
0441 int SubtractTowers::CreateNode(PHCompositeNode *topNode)
0442 {
0443   PHNodeIterator iter(topNode);
0444 
0445   // Looking for the DST node
0446   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0447   if (!dstNode)
0448   {
0449     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0450     return Fun4AllReturnCodes::ABORTRUN;
0451   }
0452 
0453   IHTowerName = m_towerNodePrefix + "_HCALIN";
0454   TowerInfoContainer *hcal_towers = findNode::getClass<TowerInfoContainer>(topNode, IHTowerName);
0455   if (m_use_towerinfo && !hcal_towers)
0456   {
0457     std::cout << PHWHERE << "Cannot find " << IHTowerName << " for creating new tower containers. Exiting" << std::endl;
0458     exit(1);
0459   }
0460 
0461   // store the new EMCal towers
0462 
0463   PHCompositeNode *emcalNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "CEMC"));
0464   if (!emcalNode)
0465   {
0466     std::cout << PHWHERE << "EMCal Node note found, doing nothing." << std::endl;
0467   }
0468   if (m_use_towerinfo)
0469   {
0470     EMTowerName = m_towerNodePrefix + "_CEMC_RETOWER_SUB1";
0471     TowerInfoContainer *test_emcal_tower = findNode::getClass<TowerInfoContainer>(topNode, EMTowerName);
0472     if (!test_emcal_tower)
0473     {
0474       if (Verbosity() > 0)
0475       {
0476         std::cout << "SubtractTowers::CreateNode : creating " << EMTowerName << " node " << std::endl;
0477       }
0478 
0479       TowerInfoContainer *emcal_towers = dynamic_cast<TowerInfoContainer *>(hcal_towers->CloneMe());
0480       PHIODataNode<PHObject> *emcalTowerNode = new PHIODataNode<PHObject>(emcal_towers, EMTowerName, "PHObject");
0481       emcalNode->addNode(emcalTowerNode);
0482     }
0483     else
0484     {
0485       std::cout << "SubtractTowers::CreateNode : " << EMTowerName << " already exists! " << std::endl;
0486     }
0487   }
0488   else
0489   {
0490     RawTowerContainer *test_emcal_tower = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER_SUB1");
0491     if (!test_emcal_tower)
0492     {
0493       if (Verbosity() > 0)
0494       {
0495         std::cout << "SubtractTowers::CreateNode : creating TOWER_CALIB_CEMC_RETOWER_SUB1 node " << std::endl;
0496       }
0497 
0498       RawTowerContainer *emcal_towers = new RawTowerContainer(RawTowerDefs::CalorimeterId::HCALIN);
0499       PHIODataNode<PHObject> *emcalTowerNode = new PHIODataNode<PHObject>(emcal_towers, "TOWER_CALIB_CEMC_RETOWER_SUB1", "PHObject");
0500       emcalNode->addNode(emcalTowerNode);
0501     }
0502     else
0503     {
0504       std::cout << "SubtractTowers::CreateNode : TOWER_CALIB_CEMC_RETOWER_SUB1 already exists! " << std::endl;
0505     }
0506   }
0507 
0508   // store the new IHCal towers
0509   PHCompositeNode *ihcalNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "HCALIN"));
0510   if (!ihcalNode)
0511   {
0512     std::cout << PHWHERE << "IHCal Node note found, doing nothing." << std::endl;
0513   }
0514   if (m_use_towerinfo)
0515   {
0516     IHTowerName = m_towerNodePrefix + "_HCALIN_SUB1";
0517     TowerInfoContainer *test_ihcal_tower = findNode::getClass<TowerInfoContainer>(topNode, IHTowerName);
0518     if (!test_ihcal_tower)
0519     {
0520       if (Verbosity() > 0)
0521       {
0522         std::cout << "SubtractTowers::CreateNode : creating " << IHTowerName << " node " << std::endl;
0523       }
0524 
0525       TowerInfoContainer *ihcal_towers = dynamic_cast<TowerInfoContainer *>(hcal_towers->CloneMe());
0526       PHIODataNode<PHObject> *ihcalTowerNode = new PHIODataNode<PHObject>(ihcal_towers, IHTowerName, "PHObject");
0527       ihcalNode->addNode(ihcalTowerNode);
0528     }
0529     else
0530     {
0531       std::cout << "SubtractTowers::CreateNode : " << IHTowerName << " already exists! " << std::endl;
0532     }
0533   }
0534   else
0535   {
0536     RawTowerContainer *test_ihcal_tower = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN_SUB1");
0537     if (!test_ihcal_tower)
0538     {
0539       if (Verbosity() > 0)
0540       {
0541         std::cout << "SubtractTowers::CreateNode : creating TOWER_CALIB_HCALIN_SUB1 node " << std::endl;
0542       }
0543 
0544       RawTowerContainer *ihcal_towers = new RawTowerContainer(RawTowerDefs::CalorimeterId::HCALIN);
0545       PHIODataNode<PHObject> *ihcalTowerNode = new PHIODataNode<PHObject>(ihcal_towers, "TOWER_CALIB_HCALIN_SUB1", "PHObject");
0546       ihcalNode->addNode(ihcalTowerNode);
0547     }
0548     else
0549     {
0550       std::cout << "SubtractTowers::CreateNode : TOWER_CALIB_HCALIN_SUB1 already exists! " << std::endl;
0551     }
0552   }
0553 
0554   // store the new OHCal towers
0555   PHCompositeNode *ohcalNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "HCALOUT"));
0556   if (!ohcalNode)
0557   {
0558     std::cout << PHWHERE << "OHCal Node note found, doing nothing." << std::endl;
0559   }
0560   if (m_use_towerinfo)
0561   {
0562     OHTowerName = m_towerNodePrefix + "_HCALOUT_SUB1";
0563     TowerInfoContainer *test_ohcal_tower = findNode::getClass<TowerInfoContainer>(topNode, OHTowerName);
0564     if (!test_ohcal_tower)
0565     {
0566       if (Verbosity() > 0)
0567       {
0568         std::cout << "SubtractTowers::CreateNode : creating " << OHTowerName << " node " << std::endl;
0569       }
0570 
0571       TowerInfoContainer *ohcal_towers = dynamic_cast<TowerInfoContainer *>(hcal_towers->CloneMe());
0572       PHIODataNode<PHObject> *ohcalTowerNode = new PHIODataNode<PHObject>(ohcal_towers, OHTowerName, "PHObject");
0573       ohcalNode->addNode(ohcalTowerNode);
0574     }
0575     else
0576     {
0577       std::cout << "SubtractTowers::CreateNode : " << OHTowerName << " already exists! " << std::endl;
0578     }
0579   }
0580   else
0581   {
0582     RawTowerContainer *test_ohcal_tower = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT_SUB1");
0583     if (!test_ohcal_tower)
0584     {
0585       if (Verbosity() > 0)
0586       {
0587         std::cout << "SubtractTowers::CreateNode : creating TOWER_CALIB_HCALOUT_SUB1 node " << std::endl;
0588       }
0589 
0590       RawTowerContainer *ohcal_towers = new RawTowerContainer(RawTowerDefs::CalorimeterId::HCALOUT);
0591       PHIODataNode<PHObject> *ohcalTowerNode = new PHIODataNode<PHObject>(ohcal_towers, "TOWER_CALIB_HCALOUT_SUB1", "PHObject");
0592       ohcalNode->addNode(ohcalTowerNode);
0593     }
0594     else
0595     {
0596       std::cout << "SubtractTowers::CreateNode : TOWER_CALIB_HCALOUT_SUB1 already exists! " << std::endl;
0597     }
0598   }
0599   return Fun4AllReturnCodes::EVENT_OK;
0600 }