packages/jetbackground/SubtractTowersCS.cc

0001 #include "SubtractTowersCS.h"
0002
0003 #include "TowerBackground.h"
0004
0005 #include <calobase/RawTower.h>
0006 #include <calobase/RawTowerContainer.h>
0007 #include <calobase/RawTowerDefs.h>
0008 #include <calobase/RawTowerGeom.h>
0009 #include <calobase/RawTowerGeomContainer.h>
0010 #include <calobase/RawTowerv1.h>
0011
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013 #include <fun4all/SubsysReco.h>
0014
0015 #include <phool/PHCompositeNode.h>
0016 #include <phool/PHIODataNode.h>
0017 #include <phool/PHNode.h>
0018 #include <phool/PHNodeIterator.h>
0019 #include <phool/PHObject.h>
0020 #include <phool/getClass.h>
0021 #include <phool/phool.h>
0022
0023 #include <fastjet/PseudoJet.hh>
0024 #include <fastjet/contrib/ConstituentSubtractor.hh>
0025
0026 // standard includes
0027 #include <cmath>
0028 #include <iostream>
0029 #include <map>
0030 #include <utility>
0031 #include <vector>
0032
0033 SubtractTowersCS::SubtractTowersCS(const std::string &name)
0034   : SubsysReco(name)
0035 {
0036 }
0037
0038 int SubtractTowersCS::InitRun(PHCompositeNode *topNode)
0039 {
0040   CreateNode(topNode);
0041
0042   return Fun4AllReturnCodes::EVENT_OK;
0043 }
0044
0045 int SubtractTowersCS::process_event(PHCompositeNode *topNode)
0046 {
0047   if (Verbosity() > 0)
0048   {
0049     std::cout << "SubtractTowersCS::process_event: entering, with _use_flow_modulation = " << _use_flow_modulation << std::endl;
0050   }
0051
0052   // pull out the tower containers and geometry objects at the start
0053   RawTowerContainer *towersEM3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER");
0054   RawTowerContainer *towersIH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
0055   RawTowerContainer *towersOH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
0056
0057   RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0058   RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0059
0060   if (Verbosity() > 0)
0061   {
0062     std::cout << "SubtractTowersCS::process_event: " << towersEM3->size() << " TOWER_CALIB_CEMC_RETOWER towers" << std::endl;
0063     std::cout << "SubtractTowersCS::process_event: " << towersIH3->size() << " TOWER_CALIB_HCALIN towers" << std::endl;
0064     std::cout << "SubtractTowersCS::process_event: " << towersOH3->size() << " TOWER_CALIB_HCALOUT towers" << std::endl;
0065   }
0066
0067   // these should have already been created during InitRun()
0068   // note that these are the CS-versions
0069   RawTowerContainer *emcal_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER_SUB1CS");
0070   RawTowerContainer *ihcal_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN_SUB1CS");
0071   RawTowerContainer *ohcal_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT_SUB1CS");
0072
0073   if (Verbosity() > 0)
0074   {
0075     std::cout << "SubtractTowersCS::process_event: starting with " << emcal_towers->size() << " TOWER_CALIB_CEMC_RETOWER_SUB1CS towers" << std::endl;
0076     std::cout << "SubtractTowersCS::process_event: starting with " << ihcal_towers->size() << " TOWER_CALIB_HCALIN_SUB1CS towers" << std::endl;
0077     std::cout << "SubtractTowersCS::process_event: starting with " << ohcal_towers->size() << " TOWER_CALIB_HCALOUT_SUB1CS towers" << std::endl;
0078   }
0079
0080   TowerBackground *towerbackground = findNode::getClass<TowerBackground>(topNode, "TowerBackground_Sub2");
0081
0082   // read these in to use, even if we don't use flow modulation in the subtraction
0083   float background_v2 = towerbackground->get_v2();
0084   float background_Psi2 = towerbackground->get_Psi2();
0085
0086   // set up constituent subtraction
0087   fastjet::contrib::ConstituentSubtractor subtractor;
0088
0089   // free parameter for the type of distance between particle i and
0090   // ghost k. There are two options: "deltaR" or "angle" which are
0091   // defined as deltaR=sqrt((y_i-y_k)^2+(phi_i-phi_k)^2) or Euclidean
0092   // angle between the momenta
0093   subtractor.set_distance_type(fastjet::contrib::ConstituentSubtractor::deltaR);
0094
0095   // free parameter for the maximal allowed distance between particle i and ghost k
0096   subtractor.set_max_distance(_DeltaRmax);
0097
0098   // free parameter for the distance measure (the exponent of particle
0099   // pt). The larger the parameter alpha, the more are favoured the
0100   // lower pt particles in the subtraction process
0101   subtractor.set_alpha(_alpha);
0102
0103   // free parameter for the density of ghosts. The smaller, the better
0104   // - but also the computation is slower.
0105   subtractor.set_ghost_area(0.01);
0106
0107   // EM layer first
0108   std::vector<fastjet::PseudoJet> backgroundProxies_EM;
0109   std::vector<fastjet::PseudoJet> fullEvent_EM;
0110   std::vector<fastjet::PseudoJet> *backgroundProxies_EM_remaining = new std::vector<fastjet::PseudoJet>;
0111
0112   // create PseudoJet version of EM towers in unsubtracted event
0113   RawTowerContainer::ConstRange begin_end_EM = towersEM3->getTowers();
0114
0115   for (RawTowerContainer::ConstIterator rtiter = begin_end_EM.first; rtiter != begin_end_EM.second; ++rtiter)
0116   {
0117     RawTower *tower = rtiter->second;
0118
0119     double this_eta = geomIH->get_tower_geometry(tower->get_key())->get_eta();
0120     double this_phi = geomIH->get_tower_geometry(tower->get_key())->get_phi();
0121     double this_E = tower->get_energy();
0122
0123     double this_pz = this_E * tanh(this_eta);
0124     double this_pt = sqrt(pow(this_E, 2) - pow(this_pz, 2));
0125     double this_px = this_pt * cos(this_phi);
0126     double this_py = this_pt * sin(this_phi);
0127
0128     fastjet::PseudoJet this_pj(this_px, this_py, this_pz, this_E);
0129
0130     fullEvent_EM.push_back(this_pj);
0131   }
0132
0133   // create all new towers
0134   for (int eta = 0; eta < geomIH->get_etabins(); eta++)
0135   {
0136     for (int phi = 0; phi < geomIH->get_phibins(); phi++)
0137     {
0138       RawTower *new_tower = new RawTowerv1();
0139       new_tower->set_energy(0);
0140       emcal_towers->AddTower(eta, phi, new_tower);
0141     }
0142   }
0143
0144   // create PseudoJet version of estimated background in EM layer
0145   for (RawTowerContainer::ConstIterator rtiter = emcal_towers->getTowers().first; rtiter != emcal_towers->getTowers().second; ++rtiter)
0146   {
0147     RawTower *tower = rtiter->second;
0148
0149     double this_eta = geomIH->get_tower_geometry(tower->get_key())->get_eta();
0150     double this_phi = geomIH->get_tower_geometry(tower->get_key())->get_phi();
0151
0152     double UE = towerbackground->get_UE(0).at(tower->get_bineta());
0153     if (_use_flow_modulation)
0154     {
0155       UE = UE * (1 + 2 * background_v2 * cos(2 * (this_phi - background_Psi2)));
0156     }
0157
0158     double this_pz = UE * tanh(this_eta);
0159     double this_pt = sqrt(pow(UE, 2) - pow(this_pz, 2));
0160     double this_px = this_pt * cos(this_phi);
0161     double this_py = this_pt * sin(this_phi);
0162
0163     fastjet::PseudoJet this_pj(this_px, this_py, this_pz, UE);
0164
0165     backgroundProxies_EM.push_back(this_pj);
0166
0167     if (Verbosity() > 5)
0168     {
0169       std::cout << " SubtractTowersCS::process_event : background tower EM estimate for eta / phi = " << tower->get_bineta() << " / " << tower->get_binphi() << ", UE = " << UE << std::endl;
0170     }
0171   }
0172
0173   // constituent subtraction
0174   std::vector<fastjet::PseudoJet> correctedEvent_EM = subtractor.do_subtraction(fullEvent_EM, backgroundProxies_EM, backgroundProxies_EM_remaining);
0175
0176   if (Verbosity() > 0)
0177   {
0178     std::cout << " SubtractTowersCS::process_event : vector lengths fullEvent_EM = " << fullEvent_EM.size() << " , backgroundProxies_EM = " << backgroundProxies_EM.size() << " , correctedEvent_EM = " << correctedEvent_EM.size() << ", backgroundProxies_EM_remaining = " << backgroundProxies_EM_remaining->size() << std::endl;
0179
0180     double E_0 = 0;
0181     double E_1 = 0;
0182     double E_2 = 0;
0183     double E_3 = 0;
0184
0185     for (auto &n : fullEvent_EM)
0186     {
0187       E_0 += n.E();
0188     }
0189     for (auto &n : backgroundProxies_EM)
0190     {
0191       E_1 += n.E();
0192     }
0193     for (auto &n : correctedEvent_EM)
0194     {
0195       E_2 += n.E();
0196     }
0197     for (auto &n : *backgroundProxies_EM_remaining)
0198     {
0199       E_3 += n.E();
0200     }
0201
0202     std::cout << "SubtractTowersCS::process_event EM : full event E - background E = " << E_0 << " - " << E_1 << " = " << E_0 - E_1 << ", subtracted E - remaining bkg E = " << E_2 << " - " << E_3 << " = " << E_2 - E_3 << std::endl;
0203   }
0204
0205   // load subtracted towers into grid
0206
0207   for (auto &n : correctedEvent_EM)
0208   {
0209     float this_eta = n.eta();
0210     float this_phi = n.phi();
0211     float this_E = n.E();
0212
0213     // look up tower by eta / phi...
0214
0215     int this_etabin = geomIH->get_etabin(this_eta);
0216     int this_phibin = geomIH->get_phibin(this_phi);
0217     RawTower *tower = emcal_towers->getTower(this_etabin, this_phibin);
0218     tower->set_energy(this_E);
0219
0220     if (Verbosity() > 5 || this_E < 0)
0221     {
0222       std::cout << " SubtractTowersCS::process_event : creating subtracted EM tower for eta / phi = " << this_eta << " / " << this_phi << " ( " << tower->get_binphi() << " , " << this_phibin << " ) , sub. E  = " << this_E << (this_E < 0 ? " -- WARNING: negative E" : "") << std::endl;
0223     }
0224   }
0225
0226   // IHCal layer
0227   std::vector<fastjet::PseudoJet> backgroundProxies_IH;
0228   std::vector<fastjet::PseudoJet> fullEvent_IH;
0229   std::vector<fastjet::PseudoJet> *backgroundProxies_IH_remaining = new std::vector<fastjet::PseudoJet>;
0230
0231   // create PseudoJet version of IH towers in unsubtracted event
0232   RawTowerContainer::ConstRange begin_end_IH = towersIH3->getTowers();
0233
0234   for (RawTowerContainer::ConstIterator rtiter = begin_end_IH.first; rtiter != begin_end_IH.second; ++rtiter)
0235   {
0236     RawTower *tower = rtiter->second;
0237
0238     double this_eta = geomIH->get_tower_geometry(tower->get_key())->get_eta();
0239     double this_phi = geomIH->get_tower_geometry(tower->get_key())->get_phi();
0240     double this_E = tower->get_energy();
0241
0242     double this_pz = this_E * tanh(this_eta);
0243     double this_pt = sqrt(pow(this_E, 2) - pow(this_pz, 2));
0244     double this_px = this_pt * cos(this_phi);
0245     double this_py = this_pt * sin(this_phi);
0246
0247     fastjet::PseudoJet this_pj(this_px, this_py, this_pz, this_E);
0248
0249     fullEvent_IH.push_back(this_pj);
0250   }
0251
0252   // create all new towers
0253   for (int eta = 0; eta < geomIH->get_etabins(); eta++)
0254   {
0255     for (int phi = 0; phi < geomIH->get_phibins(); phi++)
0256     {
0257       RawTower *new_tower = new RawTowerv1();
0258       new_tower->set_energy(0);
0259       ihcal_towers->AddTower(eta, phi, new_tower);
0260     }
0261   }
0262
0263   // create PseudoJet version of estimated background in IH layer
0264   for (RawTowerContainer::ConstIterator rtiter = ihcal_towers->getTowers().first; rtiter != ihcal_towers->getTowers().second; ++rtiter)
0265   {
0266     RawTower *tower = rtiter->second;
0267
0268     double this_eta = geomIH->get_tower_geometry(tower->get_key())->get_eta();
0269     double this_phi = geomIH->get_tower_geometry(tower->get_key())->get_phi();
0270
0271     double UE = towerbackground->get_UE(1).at(tower->get_bineta());
0272     if (_use_flow_modulation)
0273     {
0274       UE = UE * (1 + 2 * background_v2 * cos(2 * (this_phi - background_Psi2)));
0275     }
0276
0277     double this_pz = UE * tanh(this_eta);
0278     double this_pt = sqrt(pow(UE, 2) - pow(this_pz, 2));
0279     double this_px = this_pt * cos(this_phi);
0280     double this_py = this_pt * sin(this_phi);
0281
0282     fastjet::PseudoJet this_pj(this_px, this_py, this_pz, UE);
0283
0284     backgroundProxies_IH.push_back(this_pj);
0285
0286     if (Verbosity() > 5)
0287     {
0288       std::cout << " SubtractTowersCS::process_event : background tower IH estimate for eta / phi = " << tower->get_bineta() << " / " << tower->get_binphi() << ", UE = " << UE << std::endl;
0289     }
0290   }
0291
0292   // constituent subtraction
0293   std::vector<fastjet::PseudoJet> correctedEvent_IH = subtractor.do_subtraction(fullEvent_IH, backgroundProxies_IH, backgroundProxies_IH_remaining);
0294
0295   if (Verbosity() > 0)
0296   {
0297     std::cout << " SubtractTowersCS::process_event : vector lengths fullEvent_IH = " << fullEvent_IH.size() << " , backgroundProxies_IH = " << backgroundProxies_IH.size() << " , correctedEvent_IH = " << correctedEvent_IH.size() << ", backgroundProxies_IH_remaining = " << backgroundProxies_IH_remaining->size() << std::endl;
0298
0299     double E_0 = 0;
0300     double E_1 = 0;
0301     double E_2 = 0;
0302     double E_3 = 0;
0303
0304     for (auto &n : fullEvent_IH)
0305     {
0306       E_0 += n.E();
0307     }
0308     for (auto &n : backgroundProxies_IH)
0309     {
0310       E_1 += n.E();
0311     }
0312     for (auto &n : correctedEvent_IH)
0313     {
0314       E_2 += n.E();
0315     }
0316     for (auto &n : *backgroundProxies_IH_remaining)
0317     {
0318       E_3 += n.E();
0319     }
0320
0321     std::cout << "SubtractTowersCS::process_event IH : full event E - background E = " << E_0 << " - " << E_1 << " = " << E_0 - E_1 << ", subtracted E - remaining bkg E = " << E_2 << " - " << E_3 << " = " << E_2 - E_3 << std::endl;
0322   }
0323
0324   // load subtracted towers into grid
0325
0326   for (auto &n : correctedEvent_IH)
0327   {
0328     float this_eta = n.eta();
0329     float this_phi = n.phi();
0330     float this_E = n.E();
0331
0332     // look up tower by eta / phi...
0333
0334     int this_etabin = geomIH->get_etabin(this_eta);
0335     int this_phibin = geomIH->get_phibin(this_phi);
0336     RawTower *tower = ihcal_towers->getTower(this_etabin, this_phibin);
0337     tower->set_energy(this_E);
0338
0339     if (Verbosity() > 5 || this_E < 0)
0340     {
0341       std::cout << " SubtractTowersCS::process_event : creating subtracted IH tower for eta / phi = " << this_eta << " / " << this_phi << " ( " << tower->get_binphi() << " , " << this_phibin << " ) , sub. E  = " << this_E << (this_E < 0 ? " -- WARNING: negative E" : "") << std::endl;
0342     }
0343   }
0344
0345   // OHCal layer
0346   std::vector<fastjet::PseudoJet> backgroundProxies_OH;
0347   std::vector<fastjet::PseudoJet> fullEvent_OH;
0348   std::vector<fastjet::PseudoJet> *backgroundProxies_OH_remaining = new std::vector<fastjet::PseudoJet>;
0349
0350   // create PseudoJet version of OH towers in unsubtracted event
0351   RawTowerContainer::ConstRange begin_end_OH = towersOH3->getTowers();
0352
0353   for (RawTowerContainer::ConstIterator rtiter = begin_end_OH.first; rtiter != begin_end_OH.second; ++rtiter)
0354   {
0355     RawTower *tower = rtiter->second;
0356
0357     double this_eta = geomOH->get_tower_geometry(tower->get_key())->get_eta();
0358     double this_phi = geomOH->get_tower_geometry(tower->get_key())->get_phi();
0359     double this_E = tower->get_energy();
0360
0361     double this_pz = this_E * tanh(this_eta);
0362     double this_pt = sqrt(pow(this_E, 2) - pow(this_pz, 2));
0363     double this_px = this_pt * cos(this_phi);
0364     double this_py = this_pt * sin(this_phi);
0365
0366     fastjet::PseudoJet this_pj(this_px, this_py, this_pz, this_E);
0367
0368     fullEvent_OH.push_back(this_pj);
0369   }
0370
0371   // create all new towers
0372   for (int eta = 0; eta < geomOH->get_etabins(); eta++)
0373   {
0374     for (int phi = 0; phi < geomOH->get_phibins(); phi++)
0375     {
0376       RawTower *new_tower = new RawTowerv1();
0377       new_tower->set_energy(0);
0378       ohcal_towers->AddTower(eta, phi, new_tower);
0379     }
0380   }
0381
0382   // create PseudoJet version of estimated background in OH layer
0383   for (RawTowerContainer::ConstIterator rtiter = ohcal_towers->getTowers().first; rtiter != ohcal_towers->getTowers().second; ++rtiter)
0384   {
0385     RawTower *tower = rtiter->second;
0386
0387     double this_eta = geomOH->get_tower_geometry(tower->get_key())->get_eta();
0388     double this_phi = geomOH->get_tower_geometry(tower->get_key())->get_phi();
0389
0390     double UE = towerbackground->get_UE(2).at(tower->get_bineta());
0391     if (_use_flow_modulation)
0392     {
0393       UE = UE * (1 + 2 * background_v2 * cos(2 * (this_phi - background_Psi2)));
0394     }
0395
0396     double this_pz = UE * tanh(this_eta);
0397     double this_pt = sqrt(pow(UE, 2) - pow(this_pz, 2));
0398     double this_px = this_pt * cos(this_phi);
0399     double this_py = this_pt * sin(this_phi);
0400
0401     fastjet::PseudoJet this_pj(this_px, this_py, this_pz, UE);
0402
0403     backgroundProxies_OH.push_back(this_pj);
0404
0405     if (Verbosity() > 5)
0406     {
0407       std::cout << " SubtractTowersCS::process_event : background tower OH estimate for eta / phi = " << tower->get_bineta() << " / " << tower->get_binphi() << ", UE = " << UE << std::endl;
0408     }
0409   }
0410
0411   // constituent subtraction
0412   std::vector<fastjet::PseudoJet> correctedEvent_OH = subtractor.do_subtraction(fullEvent_OH, backgroundProxies_OH, backgroundProxies_OH_remaining);
0413
0414   if (Verbosity() > 0)
0415   {
0416     std::cout << " SubtractTowersCS::process_event : vector lengths fullEvent_OH = " << fullEvent_OH.size() << " , backgroundProxies_OH = " << backgroundProxies_OH.size() << " , correctedEvent_OH = " << correctedEvent_OH.size() << ", backgroundProxies_OH_remaining = " << backgroundProxies_OH_remaining->size() << std::endl;
0417
0418     double E_0 = 0;
0419     double E_1 = 0;
0420     double E_2 = 0;
0421     double E_3 = 0;
0422
0423     for (auto &n : fullEvent_OH)
0424     {
0425       E_0 += n.E();
0426     }
0427     for (auto &n : backgroundProxies_OH)
0428     {
0429       E_1 += n.E();
0430     }
0431     for (auto &n : correctedEvent_OH)
0432     {
0433       E_2 += n.E();
0434     }
0435     for (auto &n : *backgroundProxies_OH_remaining)
0436     {
0437       E_3 += n.E();
0438     }
0439
0440     std::cout << "SubtractTowersCS::process_event OH : full event E - background E = " << E_0 << " - " << E_1 << " = " << E_0 - E_1 << ", subtracted E - remaining bkg E = " << E_2 << " - " << E_3 << " = " << E_2 - E_3 << std::endl;
0441   }
0442
0443   // load subtracted towers into grid
0444
0445   for (auto &n : correctedEvent_OH)
0446   {
0447     float this_eta = n.eta();
0448     float this_phi = n.phi();
0449     float this_E = n.E();
0450
0451     // look up tower by eta / phi...
0452
0453     int this_etabin = geomOH->get_etabin(this_eta);
0454     int this_phibin = geomOH->get_phibin(this_phi);
0455     RawTower *tower = ohcal_towers->getTower(this_etabin, this_phibin);
0456     tower->set_energy(this_E);
0457
0458     if (Verbosity() > 5 || this_E < 0)
0459     {
0460       std::cout << " SubtractTowersCS::process_event : creating subtracted OH tower for eta / phi = " << this_eta << " / " << this_phi << " ( " << tower->get_binphi() << " , " << this_phibin << " ) , sub. E  = " << this_E << (this_E < 0 ? " -- WARNING: negative E" : "") << std::endl;
0461     }
0462   }
0463
0464   // cleanup
0465
0466   delete backgroundProxies_EM_remaining;
0467   delete backgroundProxies_IH_remaining;
0468   delete backgroundProxies_OH_remaining;
0469
0470   if (Verbosity() > 0)
0471   {
0472     std::cout << "SubtractTowersCS::process_event: ending with " << emcal_towers->size() << " TOWER_CALIB_CEMC_RETOWER_SUB1CS towers" << std::endl;
0473
0474     float EM_old_E = 0;
0475     float EM_new_E = 0;
0476     {
0477       RawTowerContainer::ConstRange begin_end_EM_1 = towersEM3->getTowers();
0478       for (RawTowerContainer::ConstIterator rtiter = begin_end_EM_1.first; rtiter != begin_end_EM_1.second; ++rtiter)
0479       {
0480         RawTower *tower = rtiter->second;
0481         EM_old_E += tower->get_energy();
0482       }
0483     }
0484     {
0485       RawTowerContainer::ConstRange begin_end_EM_2 = emcal_towers->getTowers();
0486       for (RawTowerContainer::ConstIterator rtiter = begin_end_EM_2.first; rtiter != begin_end_EM_2.second; ++rtiter)
0487       {
0488         RawTower *tower = rtiter->second;
0489         EM_new_E += tower->get_energy();
0490       }
0491     }
0492     std::cout << "SubtractTowersCS::process_event: old / new total E in EM layer = " << EM_old_E << " / " << EM_new_E << std::endl;
0493
0494     std::cout << "SubtractTowersCS::process_event: ending with " << ihcal_towers->size() << " TOWER_CALIB_HCALIN_SUB1CS towers" << std::endl;
0495
0496     float IH_old_E = 0;
0497     float IH_new_E = 0;
0498     {
0499       RawTowerContainer::ConstRange begin_end_EM_3 = towersIH3->getTowers();
0500       for (RawTowerContainer::ConstIterator rtiter = begin_end_EM_3.first; rtiter != begin_end_EM_3.second; ++rtiter)
0501       {
0502         RawTower *tower = rtiter->second;
0503         IH_old_E += tower->get_energy();
0504       }
0505     }
0506     {
0507       RawTowerContainer::ConstRange begin_end_EM_4 = ihcal_towers->getTowers();
0508       for (RawTowerContainer::ConstIterator rtiter = begin_end_EM_4.first; rtiter != begin_end_EM_4.second; ++rtiter)
0509       {
0510         RawTower *tower = rtiter->second;
0511         IH_new_E += tower->get_energy();
0512       }
0513     }
0514     std::cout << "SubtractTowersCS::process_event: old / new total E in IH layer = " << IH_old_E << " / " << IH_new_E << std::endl;
0515
0516     std::cout << "SubtractTowersCS::process_event: ending with " << ohcal_towers->size() << " TOWER_CALIB_HCALOUT_SUB1CS towers" << std::endl;
0517
0518     float OH_old_E = 0;
0519     float OH_new_E = 0;
0520     {
0521       RawTowerContainer::ConstRange begin_end_EM_5 = towersOH3->getTowers();
0522       for (RawTowerContainer::ConstIterator rtiter = begin_end_EM_5.first; rtiter != begin_end_EM_5.second; ++rtiter)
0523       {
0524         RawTower *tower = rtiter->second;
0525         OH_old_E += tower->get_energy();
0526       }
0527     }
0528     {
0529       RawTowerContainer::ConstRange begin_end_EM_6 = ohcal_towers->getTowers();
0530       for (RawTowerContainer::ConstIterator rtiter = begin_end_EM_6.first; rtiter != begin_end_EM_6.second; ++rtiter)
0531       {
0532         RawTower *tower = rtiter->second;
0533         OH_new_E += tower->get_energy();
0534       }
0535     }
0536     std::cout << "SubtractTowersCS::process_event: old / new total E in OH layer = " << OH_old_E << " / " << OH_new_E << std::endl;
0537   }
0538
0539   if (Verbosity() > 0)
0540   {
0541     std::cout << "SubtractTowersCS::process_event: exiting" << std::endl;
0542   }
0543
0544   return Fun4AllReturnCodes::EVENT_OK;
0545 }
0546
0547 int SubtractTowersCS::CreateNode(PHCompositeNode *topNode)
0548 {
0549   PHNodeIterator iter(topNode);
0550
0551   // Looking for the DST node
0552   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0553   if (!dstNode)
0554   {
0555     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0556     return Fun4AllReturnCodes::ABORTRUN;
0557   }
0558
0559   // store the new EMCal towers
0560   PHCompositeNode *emcalNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "CEMC"));
0561   if (!emcalNode)
0562   {
0563     std::cout << PHWHERE << "EMCal Node note found, doing nothing." << std::endl;
0564   }
0565
0566   RawTowerContainer *test_emcal_tower = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER_SUB1CS");
0567   if (!test_emcal_tower)
0568   {
0569     if (Verbosity() > 0)
0570     {
0571       std::cout << "SubtractTowersCS::CreateNode : creating TOWER_CALIB_CEMC_RETOWER_SUB1CS node " << std::endl;
0572     }
0573
0574     RawTowerContainer *emcal_towers = new RawTowerContainer(RawTowerDefs::CalorimeterId::HCALIN);
0575     PHIODataNode<PHObject> *emcalTowerNode = new PHIODataNode<PHObject>(emcal_towers, "TOWER_CALIB_CEMC_RETOWER_SUB1CS", "PHObject");
0576     emcalNode->addNode(emcalTowerNode);
0577   }
0578   else
0579   {
0580     std::cout << "SubtractTowersCS::CreateNode : TOWER_CALIB_CEMC_RETOWER_SUB1CS already exists! " << std::endl;
0581   }
0582
0583   // store the new IHCal towers
0584   PHCompositeNode *ihcalNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "HCALIN"));
0585   if (!ihcalNode)
0586   {
0587     std::cout << PHWHERE << "IHCal Node note found, doing nothing." << std::endl;
0588   }
0589
0590   RawTowerContainer *test_ihcal_tower = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN_SUB1CS");
0591   if (!test_ihcal_tower)
0592   {
0593     if (Verbosity() > 0)
0594     {
0595       std::cout << "SubtractTowersCS::CreateNode : creating TOWER_CALIB_HCALIN_SUB1CS node " << std::endl;
0596     }
0597
0598     RawTowerContainer *ihcal_towers = new RawTowerContainer(RawTowerDefs::CalorimeterId::HCALIN);
0599     PHIODataNode<PHObject> *ihcalTowerNode = new PHIODataNode<PHObject>(ihcal_towers, "TOWER_CALIB_HCALIN_SUB1CS", "PHObject");
0600     ihcalNode->addNode(ihcalTowerNode);
0601   }
0602   else
0603   {
0604     std::cout << "SubtractTowersCS::CreateNode : TOWER_CALIB_HCALIN_SUB1CS already exists! " << std::endl;
0605   }
0606
0607   // store the new OHCal towers
0608   PHCompositeNode *ohcalNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "HCALOUT"));
0609   if (!ohcalNode)
0610   {
0611     std::cout << PHWHERE << "OHCal Node note found, doing nothing." << std::endl;
0612   }
0613
0614   RawTowerContainer *test_ohcal_tower = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT_SUB1CS");
0615   if (!test_ohcal_tower)
0616   {
0617     if (Verbosity() > 0)
0618     {
0619       std::cout << "SubtractTowersCS::CreateNode : creating TOWER_CALIB_HCALOUT_SUB1CS node " << std::endl;
0620     }
0621
0622     RawTowerContainer *ohcal_towers = new RawTowerContainer(RawTowerDefs::CalorimeterId::HCALOUT);
0623     PHIODataNode<PHObject> *ohcalTowerNode = new PHIODataNode<PHObject>(ohcal_towers, "TOWER_CALIB_HCALOUT_SUB1CS", "PHObject");
0624     ohcalNode->addNode(ohcalTowerNode);
0625   }
0626   else
0627   {
0628     std::cout << "SubtractTowersCS::CreateNode : TOWER_CALIB_HCALOUT_SUB1CS already exists! " << std::endl;
0629   }
0630
0631   return Fun4AllReturnCodes::EVENT_OK;
0632 }