packages/jetbackground/SubtractTowers.cc

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